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I.         INTRODUCTION 

A.      PROBLEM  STATEMENT  AND  OVERVIEW 

The  demand  to  be  able  to  intercept  communication  signals  is  increasing,  as 
interception  is  a  basic  investigation  tool  for  civilian,  intelligence  and  military  author- 
ities. The  task  of  intercepting  a  communication  signal  can  be  summarized  by  detect- 
ing the  signal's  presence  in  additive  noise,  classifying  the  modulation  type,  estimating 
the  control  parameters  required  for  reception,  decoding  the  data,  and  decrypting  the 
information  content.  Sometimes,  the  procedure  is  stopped  at  the  detection  (identifi- 
cation) step  such  as  in  geo-location  or  spectrum  monitoring.  Intelligence  applications 
aim  at  the  final  step  which  is  to  extract  the  information  content. 

Communication  signals  can  utilize  a  great  diversity  of  digital  modulation  tech- 
niques. Among  the  digital  techniques,  the  spread  spectrum  signal  (SS)  is  widely  used. 
Frequency  hopping  signals  (FH)  are  a  subset  of  SS,  and  are  used  primarily  in  this 
dissertation.  In  the  open  literature,  many  signal  processing  tools  are  available  for  the 
interception  task,  and  among  those  are  the  correlation  analysis  and  wavelet  analysis. 
Each  of  these  tools  has  been  used  independently  in  the  signal  processing  and  commu- 
nication area.  In  this  work  we  will  address  the  merging  of  the  wavelet  and  correlation 
concept  to  detect,  classify  and  estimate  signal  parameters. 

The  FH  signal  is  a  non-stationary  process.  Hence  a  specific  type  of  two- 
dimensional  correlation  function,  called  the  instantaneous  correlation  function,  is  se- 
lected, as  it  accommodates  the  time- varying  nature  of  the  signal  of  interest.  Ap- 
plication of  wavelet  analysis  to  correlation  functions  is  a  new  area  and  is  still  in 
the  exploratory  stage.  The  automation  of  the  interception  and  exploitation  of  dig- 
ital communication  signals  is  the  final  goal.  This  work  addresses  the  interception, 
identification,  classification,  and  parameter  extraction  of  FH  signals. 

We  start  by  investigating  the  wavelet  transform  of  different  types  of  the  cor- 
relation functions.   Then,  we  select  the  instantaneous  correlation  function  (ICF)  as 


the  candidate  correlation  representation.  We  derive  the  instantaneous  correlation 
function  of  the  FH  signal  and  analyze  its  wavelet  transform  obtained  by  using  the 
Morlet  wavelet.  The  wavelet  transform  of  the  instantaneous  correlation  function  is 
a  3-dimensional  (3-D)  surface.  We  partition  the  3-D  surface  into  a  number  of  2-D 
surfaces,  where  each  corresponds  to  one  wavelet  scale  (called  the  wavelet  surface). 

We  address  the  interception  problem  in  two  approaches.  First,  we  visually 
inspect  the  2-D  wavelet  surfaces  to  identify  and  classify  the  structure  of  the  FH 
signal  and  obtain  an  estimate  for  the  hop  time  interval.  Second,  we  apply  a  proposed 
processing  scheme  to  estimate  the  hop  start/stop  time,  the  hop-scale  pattern,  and  the 
hop  frequency.  The  extraction  of  the  hop  start /stop  time  is  addressed  using  an  edge 
detection  approach  by  applying  a  compass  operator,  a  technique  which  is  well  known 
in  the  image  processing  area.  The  hop-scale  pattern  is  obtained  by  applying  an  energy 
analysis.  The  energy  analysis  assigns  a  scale  index  (called  the  proper  scale)  to  each 
hop.  The  proper  scale  of  each  hop  is  that  scale  which  has  the  greatest  energy  content. 
The  sequence  of  proper  scales,  representing  the  hop  sequence,  is  called  the  hop-scale 
pattern.  The  frequency  of  each  hop  can  be  extracted  from  the  wavelet  surface  at 
the  proper  scale.  Thus,  the  identification  and  classification  of  the  FH  signal  may  be 
accomplished  as  follows: 

1.  Based  on  the  hop-scale  pattern: 

If  two  or  more  wavelet  scales  are  applied,  the  hop-scale  pattern  of  the  FH 
signal  will  be  different  from  the  hop-scale  patterns  of  other  digital  modulation 
signals. 

2.  Based  on  the  frequency  diversity: 

If  all  hop  frequencies  reside  in  one  scale  then  the  FH  signal  will  have  different 
frequency  components  as  a  function  of  the  hop  intervals. 

This  dissertation  is  organized  into  eight  chapters.  In  Chapter  I,  we  introduce 
the  problem  and  review  related  work  with  a  focus  on  the  interception  of  FH  signals. 
In  Chapter  II,  we  introduce  the  most  significant  signal  analysis  tools.  In  Chapter  III, 
we  derive  the  wavelet  response  to  correlation  functions.  In  Chapter  IV,  we  derive  the 


ICF  for  the  FH  signal.  In  Chapter  V,  we  perform  the  wavelet  analysis  of  the  ICF  of 
the  FH  signal.  In  Chapter  VI,  we  analyze  the  processing  scheme.  In  Chapter  VII, 
we  provide  the  simulation  results.  Finally,  Chapter  VIII  contains  conclusions  and 
recommendations  for  future  work. 

B.      INTERCEPTION  OF  DIGITAL  COMMUNICATION 
SIGNALS 
1.       Introduction 

Interception  of  communication  signals  is  of  interest  to  a  wide  range  of  appli- 
cations in  surveillance,  intelligence,  reconnaissance,  geo-location,  spectral  monitoring 
and  jamming  [Ref.  1]. 

Digital  communication  systems  can  use  a  large  diversity  of  modulation  tech- 
niques. Examples  of  these  techniques  are:  ASK,  BPSK,  BFSK,  QAM,  MPSK  and 
MFSK.  Interception  of  digital  communication  signals  consists  of  detection  (identifi- 
cation), classification  (feature  extraction),  parameter  estimation,  decoding,  and  de- 
crypting. Through  out  this  dissertation  we  will  use  the  following  definitions: 

1.  Detection  or  identification  is  the  process  of  discriminating  between  noise  and 
a  digitally  modulated  signal. 

2.  Classification  (feature  extraction)  is  the  process  of  discriminating  between  dif- 
ferent modulation  types,  by  obtaining  significant  parameters  or  signal  constel- 
lation which  uniquely  specify  the  modulation  type. 

3.  Parameter  estimation  is  the  process  of  extracting  the  control  parameters  of 
the  signal. 

4.  Decoding  is  the  process  of  detecting  and  correcting  the  errors  of  the  commu- 
nication channel. 

5.  Decrypting  is  the  process  of  restoring  the  original  information  content  of  the 
transmitted  messages. 


2.  Interception  Tools 

A  large  number  of  open-literature  publications  address  the  interception  of 
digital  communication  signals.  The  different  signal  processing  tools  used  for  signal 
interception  may  be  categorized  as  follows: 

1.  Basic  tools:  Spectral  analysis,  correlation  analysis,  and  parametric  modeling. 

2.  Linear  tools:  Linear  transforms  including  the  wavelet  transform. 

3.  Non-linear  tools:  Higher-order  spectra,  spectral  correlation,  and  cyclic-feature 
processing. 

4.  Other  tools:  Eigen-analysis,  singular-value  decomposition,  and  stochastic  res- 
onance. 

There  are  two  different  approaches  to  address  the  interception.  The  first  approach 
addresses  the  interception  as  a  decision-making  problem  applied  to  a  set  of  alternative 
hypotheses.  This  requires  statistical  decision  theory  and  hence  a  statistical  description 
of  the  alternatives.  The  second  approach  handles  the  interception  as  a  statistical 
pattern-recognition  problem. 

3.  Digital  Signal  Modulation 

Different  tools  can  be  used  for  intercepting  digital  modulations.  For  ASK, 
QAM,  BPSK,  QPSK,  and  FSK,  applications  are  addressed  in  references  [Ref.  2]  to 
[Ref.  14],  [Ref.  17]  and  [Ref.  32].  Since  the  FSK  modulation  is  related  to  the  FH 
signal,  the  interception  tools  for  FSK  are  briefly  introduced  as  follows: 
Spectral  correlation  and  visual  inspection  is  applied  for  identification  of  FSK  signals 
in  [Ref.  2].  Higher-order  moments  are  applied  for  FSK  signals  in  [Ref.  3].  Wavelet 
analysis  of  the  FSK  signal  is  applied  for  timing  extraction  in  [Ref.  5,  17]  where  it  is 
shown  that  the  frequency  transitions  of  the  FSK  signal  are  related  to  inflection  points 
contained  in  the  scalogram. 


4.        Spread  Spectrum  Signals 

Interception  of  spread  spectrum  signals  is  addressed  in  [Ref.  18]  to  [Ref.  35]. 
There  are  five  conventional  techniques  to  intercept  frequency  hopped  (FH)  signals 
[Ref.  18,  19].  These  are: 

1.  Wideband  energy  detectors  (wide  band  radiometer), 

2.  Optimum-multichannel  FH  pulse-matched  filters, 

3.  Maximum-channel  filter-bank  combiners, 

4.  Optimum  partial-band  FH  pulse-matched  detectors,  and 

5.  Partial-band  filter-bank  combiners. 

These  techniques  differ  mainly  in  two  points.  The  first  difference  is  in  the  bandwidth 
of  the  interception  filter(s)  relative  to  the  bandwidth  of  the  FH  signal.  The  second 
difference  is  in  the  number  of  parallel  channels  relative  to  the  number  of  hopping 
frequencies.  Consequently,  they  differ  in  the  required  minimum  signal-to-noise  ratio 
(SNR)  for  acceptable  performance. 

Time  or  frequency  uncertainty  can  be  minimized  by  using  overlapping  tech- 
niques [Ref.  20].  A  likelihood-ratio-test  detector  and  a  channelized  receiver  are  widely 
accepted  as  the  optimum  system  for  detecting  FH  signals  [Ref.  20]  to  [Ref.  24].  The 
maximum-likelihood  detector  (ML),  using  a  bank  of  correlators,  is  addressed  in  [Ref. 
22].  It  assumes  prior  knowledge  of  the  hopping  times  and  frequencies,  and  one  corre- 
lator is  designated  to  each  primary  frequency  region.  Results  show  that  the  coherent 
ML  scheme  gives  a  probability  of  detection  of  0.5  at  an  SNR  of  4.5  dB  for  a  probabil- 
ity of  false  alarm  of  10-9.  For  the  noncoherent  ML  scheme,  this  performance  requires 
an  SNR  of  about  5.9  dB. 

A  generalized  likelihood  ratio  test  for  a  multiple  hop  observation  interval  is 
addressed  in  [Ref.  23].  It  assumes  prior  knowledge  of  hop  frequencies,  hop  rate, 
and  hop  times.  Using  an  observation  interval  of  103  hops  and  a  probability  of  false 
alarm  of  10~3,  a  probability  of  miss  of  10~5  is  achieved  at  an  SNR  of  8  dB.  The  same 
performance  requires  an  SNR  of  5  dB  using  an  observation  interval  of  105  hops. 


Applying  wavelet  analysis  to  the  interception  of  FH  signals  is  addressed  in 
[Ref.  35].  FH  classification  is  based  on  locating  transition  spikes  due  to  frequency 
transitions.  The  authors  suggest  using  the  STFT  for  the  hop  frequency  estimation 
instead  of  the  wavelet  transform  because  wavelets  have  varying  frequency  resolution 
and  the  FH  signal  has  a  wide  bandwidth. 


II.         FOURIER  ANALYSIS  AND  WAVELETS 

Signal  analysis  (expansion,  decomposition,  or  transformation)  is  a  method 
used  to  represent  time  signals  as  a  linear  combination  of  elementary  building  blocks 
(or  elementary  basis  functions).  This  representation  is  vital  to  the  area  of  signal 
processing.  Each  expansion  is  defined  by  its  basis  functions,  or  equivalently,  its 
basic  building  blocks.  There  are  numerous  linear  expansion  methods.  Well-known 
examples  are  the  Shannon  expansion,  the  Karhunen-Loeve  expansion,  the  Gram- 
Schmidt  expansion,  the  eigen-analysis,  and  the  most  popular,  the  Fourier  analysis. 
In  this  chapter,  we  summarize  the  evolution  of  the  analysis  tools  in  three  steps.  First, 
we  introduce  the  Fourier  analysis  as  used  for  stationary  signals.  Then,  we  introduce 
time-frequency  distributions,  which  can  be  used  to  represent  time-varying  signals. 
Finally,  we  introduce  the  concept  of  the  wavelet  analysis  to  represent  time-varying 
signals  over  a  time-scale  plane. 

A.      FOURIER  ANALYSIS 

The  major  reference  for  this  section  is  [Ref.  36].  The  Fourier  transform  (FT), 
also  called  Fourier  analysis,  is  the  most  popular  signal  decomposition  (expansion). 
The  Fourier  transform  is  used  to  represent  a  stationary  process  by  decomposing  it  into 
sinusoidal  or  complex  exponential  components.  Here,  a  stationary  process  is  defined 
as  a  process  whose  spectrum  is  not  varying  with  time.  Hence,  Fourier  techniques  work 
well  and  allow  successful  frequency  localization  of  the  spectral  components.  When 
a  non-stationary  process  is  present,  its  time-varying  spectrum  requires  an  additional 
dimension  (i.e.,  time).  The  major  problem  with  the  classical  Fourier  analysis  is  that 
it  uses  an  infinitely  long  sinusoidal  or  complex  exponential  basis  of  functions.  This 
makes  time  localization  impossible. 

By  classical  Fourier  analysis  we  mean  the  Fourier  series  (FS),  the  continuous- 
time  Fourier  transform  (FT),  the  discrete-time  Fourier  series  (DFS),  and  the  discrete- 


time  Fourier  transform  (DTFT). 

The  discrete  Fourier  transform  (DFT),  or  its  fast  implementation,  the  fast 
Fourier  transform  (FFT),  uses  finite  integration  time.  The  Fourier  transform  cannot 
give  a  time  resolution  better  than  the  integration  interval.  To  cope  with  the  re- 
quirement of  tracking  time  evolution  and  to  provide  time  localization,  the  short-time 
Fourier  transform  (STFT)  was  introduced.  The  STFT  uses  a  window  function  which 
affects  the  time  resolution;  that  is,  the  longer  the  window  interval  the  worse  is  the 
time  resolution. 

1.       Fourier  Series 

Given  a  periodic  continuous  time  signal  x(t),  with  period  T,  the  Fourier  series 
(FS)  expansion  is  given  by  : 

00 

*(*)=    E   CWe*2***',  (II.l) 

k=— oo 

where 

C(k)  =  i  /  x{t)e-j2*foktdt  (II.2) 

T  Jt 

and  f0  =  £. 

Note  that  the  expansion  coefficients  C{k)  are  evaluated  at  integer  values  of  k. 
Parseval's  relation  for  this  power  signal  (  recall  that  power  signals  have  finite  average 
power  and  infinite  energy)  is  given  by  : 

P*  =  fj  \m2dt=    £    |C(fc)|2,  (II.3) 

k=-oo 

where  Px  is  the  average  power  of  the  signal.  The  signal  has  to  satisfy  the  Dirichlet 
conditions  to  have  a  valid  Fourier  series  expansion.  These  conditions  can  be  summa- 
rized as:  the  signal  has  a  finite  number  of  discontinuities,  a  finite  number  of  maxima 
and  minima,  and  it  is  absolutely  integrable  (over  one  period). 
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2.       Fourier  Transform 

A  non-periodic  continuous  time  signal,  x(t),  can  be  represented  as 

/oo 
xuy2*ftdf,  (n.4) 

-00 

where 

/oo 
x{t)e-'2*ftdt.  (II.5) 

-00 

Note  that  both  the  signal  and  the  transform  are  continuous  functions  of  time  and 
frequency,  respectively. 

The  Dirichlet  conditions  have  to  be  met  by  the  non-periodic  signal  with  one 
modification;  that  is,  the  time  support  of  the  signal  is  considered  instead  of  the  period 
in  the  above-addressed  Dirichlet  conditions. 

Parseval's  relation  is  given  by 

/OO  /-OO 

\x(t)\2dt  =         \X(f)\2df  (II.6) 

-00  •/— 00 

where  Ex  is  the  total  energy  of  the  signal. 

Short-Time  Fourier  Transform.  Because  of  the  lack  of  time  localization,  non- 
stationary  signals  have  time-varying  spectra  which  cannot  be  represented  by  the  clas- 
sical Fourier  transform.  The  time  evolution  of  the  time- varying  spectra  needs  to  be 
considered.  Hence,  there  is  a  need  for  an  expansion  in  time  and  frequency  such  as  the 
short-time  Fourier  transform  (STFT).  This  transform  windows  the  signal  around  a 
certain  time  instant,  performs  the  frequency  domain  analysis,  and  repeats  the  process 
at  other  time  instants.  Note  that  it  is  assumed  that  the  windowed  signal  will  have 
a  non-time-varying  spectrum  (local  stationarity)  within  the  time  window.  Therefore, 
the  STFT  for  a  continuous  signal  x(t)  is  given  by 

/oo 
x{t)g*{t-l)e-JlJtdt,  (II.7) 

-oo 

where  g(t)  is  the  window  function.  X(u,  I)  is  the  spectral  description  of  x(t)  when  the 
window  is  centered  at  the  time  I.  When  the  window  g(t)  has  a  Gaussian  shape,  the 


STFT  is  called  the  Gabor  transform.  The  STFT  depends  on  a  window  function  with 
a  fixed  width  (in  time).  Recall  that  the  time  localization  and  frequency  localization 
of  the  STFT  are  controlled  by  the  window  width.  Consequently,  the  STFT  offers 
a  fixed-time  resolution  and  a  related  fixed-frequency  resolution,  which  results  in  a 
uniform  tiling  in  the  time-frequency  plane. 

The  time  and  frequency  resolution  are  governed  by  the  uncertainty  principle, 
which  will  be  addressed  shortly.  For  details,  see  [Ref.  38].  The  STFT  for  discrete 
time  signals  is  defined  as 

oo 

X{k,  n)  =  £  g{m  -  n)x{m)e-j2*kM ,  (II.8) 

m=oo 

where  k  =  0,1, ...,M -  1. 

3.  Discrete-Time  Fourier  Series 

The  discrete-time  Fourier  series  (DTFS)  for  a  periodic  discrete  signal  x(n) 
with  period  N  is  given  by 

x(n)  =  £  C(k)ej2nk%,  (11.9) 

fc=0 

and 

CM^E^Me-'2***,  (11.10) 

iV    71=0 

where  C(k)  represents  the  expansion  coefficient  at  the  discrete  frequency  Wk  =  2-Kk^. 
Parseval's  relation  is  given  by 

P.-EIOWP— sEW")!1  (nil) 

fc=0  iy   n=0 

where  Px  is  the  average  power  of  the  infinite  energy  periodic  signal  x(n).  The  Dirichlet 
conditions  must  be  met  by  x(n). 

4.  Discrete-Time  Fourier  Transform 

The  discrete-time  Fourier  transform  (DTFT)  for  finite  length  non-periodic 
discrete  signals  is  given  by 


x 


{n)  =  h  Lx{u))ejuinduj'  (IL12) 
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where 

oo 

X(u)  =    £   x(n)e~jwn-  (11.13) 

k=—oo 

Note  that  the  discrete  frequency  uj  ranges  from  —  ir  to  tt.  Parseval's  relation 

is  given  by 

oo  -1      r 

Ex=    £    \x(n)\2  =  -     \X(co)\2<k,:  (11.14) 

n=_oo  2tt  J 

where  Ex  is  the  total  energy  of  the  signal. 

Discrete  Fourier  Transform.  A  finite-length  non-periodic  discrete  signal  has  a 
continuous-frequency-domain  representation  X(u).  For  digital  signal  processing  it  is 
easier  to  represent  X{u)  by  its  discrete-frequency  samples.  This  leads  to  the  discrete 
Fourier  transform  (DFT).  The  DFT  is  given  by 

x(n)^^-j:X(k)W^  (11.15) 

7V    fc=0 

and 

X(k)  =J*  x(n)Wfr,  (11.16) 

n=0 

where  Wn  =  e~JN2L  and  n,  k  =  0, 1,  •  •  • ,  N  —  1.  The  DFT  requires  N2  complex 
multiplication  operations  and  N(N—1)  complex  addition  operations.  The  fast  Fourier 
transform  (FFT)  implements  the  DFT  with  fewer  multiplications.  The  FFT  has 
computation  complexity  (number  of  multiplications)  y/o^A7"  if  N  is  power  of  2. 

B.      TIME-FREQUENCY  DISTRIBUTIONS 
1.       Introduction 

A  time-frequency  distribution  (TFD)  describes  non-stationary  signals  by  dis- 
playing the  energy  density  over  time  and  frequency  simultaneously.  Let  the  signal 
and  its  Fourier  transform  be  denoted  by  s(t)  and  S(uj),  respectively.  The  following 
terms  and  definitions  [Ref.  38]  are  useful  in  TFD  discussions: 

•  Energy  density  or  instantaneous  power  is  given  by 

\s(t)\2.  (11.17) 
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•  Total  energy  of  the  signal  is  given  by 

f\s(t)\2dt  (11.18) 

•  Mean  time  is  given  by 

(t)=  ft\s{t)\2dt.  (11.19) 

•  Time  variance  is  given  by 

o\  =  /  (t  -  (t))2  \s(t)\2  dt,  (11.20) 

where  of  is  an  indication  of  the  duration  of  the  signal. 

•  Energy  density  spectrum  is  given  by 

\S(u)\2.  (11.21) 

It  is  the  intensity  per  Hertz  near  the  frequency  to. 

•  Parseval's  theorem  states  that 

f  \S(u)\2  duj  =  f  \s{t)\2dt.  (11.22) 

•  Mean  frequency  is  given  by 

(u)=  fcj\S(u;)\2du.  (11.23) 

•  Frequency  variance  is  given  by 

a^jiuj-i^flSiu)]2^,  (11.24) 

where  o£  is  an  indication  of  the  rms  bandwidth. 

•  The  uncertainty  principle: 

It  is  a  fundamental  statement  regarding  the  Fourier  transform  pairs.  It  governs 
the  relationship  between  the  spread  of  any  signal  in  the  time  domain  and  in 
the  frequency  domain.  This  relationship  is  given  by  Oto^  >  1/2.  This  means 
that  there  is  a  trade-off  time  localization  and  frequency  localization. 
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and 


2.       Fundamental  Properties  of  Time-Frequency 
Distributions 

The  major  reference  for  these  properties  is  [Ref.  38].  Let  p(t:ui)  be  the  inten- 
sity of  the  signal  s(t)  at  time  t  and  frequency  ui.  Then  p(t,u)At/±.uj  is  the  fractional 
energy  in  the  cell  AtAcu  at  time  t  and  frequency  cu.  Basic  properties  of  the  TFD  are 
given  by: 

•  Marginal  property:  Summing  up  the  energy  distribution  for  all  frequencies  at 
a  particular  time  gives  the  instantaneous  energy.  Summing  over  all  times  at  a 
particular  frequency  gives  the  energy  density  spectrum,  i.e., 

[p{t,to)dw  =  \s(t)\2  (11.25) 

[p(t,u)dt=  \S(uj)\2.  (11.26) 

Note:    Any  joint  distribution  that  yields  the  correct  marginal  property  is  con- 
sistent with  the  uncertainty  principle. 

•  Total  energy  property:  The  total  energy  of  the  distribution  should  be  equal 
to  the  total  energy  of  the  signal. 

E  =  f  [p{t:Lo)dujdt  =  f  \s(t)\2dt  =  f  \S(cu)\2(Lj  (11.27) 

Note:    If  the  distribution  satisfies  the  marginal  property  then  the  total  energy 
property  is  automatically  satisfied. 

•  Time-shift-invariance:  The  shifted  signal  s(t  —  t0)  will  have  the  distribution 

p(t-  t0,u>). 

•  Frequency-shift-invariance:  The  modulated  signal  s(/)e^  will  have  p(t,cj  —  £) 
as  its  distribution. 

•  Time  and  frequency  shift  in  variance:  The  signal  s(t  —  to)eJ^  will  have  p(t  — 
to^  —  0  as  its  distribution. 

•  Linear  scaling  property:  The  signal  ssc(t)  =  y/as(at)  will  have  the  distribution 
psc(t,  to)  =  p(oi,  a>/a),  which  will  satisfy  the  marginal  property  of  the  scaled 
signal  ssc(t). 

•  Weak  finite-support  property:  If  s(t)  has  zero  value  outside  the  time  interval 
(£i,  £2)  or  outside  the  frequency  range  (a^u^),  so  does  the  distribution  p(t,u). 


13 


•  Strong  finite-support  property:  If  s(t)  has  zero  value  at  any  time  to  or  at  any 
frequency  a>o,  then  p(t,u)  will  be  zero  for  any  time  to  or  frequency  u0. 

Note:  It  is  impossible  for  the  signal  f(t)  to  be  both  time  and  frequency 
limited.  Therefore,  if  a  distribution  satisfies  the  weak  finite-support  property, 
it  cannot  be  limited  to  a  finite  region  of  time-frequency  plane. 

3.       Wigner-Ville  Distribution 

The  major  reference  of  this  topic  is  [Ref.  38,  39].  The  reason  the  Wigner-Ville 
distribution  (WVD)  is  discussed  in  detail  is  that  we  relate  the  WVD  and  the  wavelet 
transform  of  the  instantaneous  correlation  function.  Another  reason  is  that  the  WVD 
is  typically  used  to  represent  a  non-stationary  process. 

The  WVD,  in  terms  of  the  signal  s(t)  or  its  spectrum  5(u>),  is  defined  by: 


and 


W(t,w)  = 


iir^m^w    <w> 


4.       Properties  of  the  Wigner-Ville  Distribution 

The  WVD  has  the  following  properties: 

•  The  WVD  is  always  real,  even  if  s(t)  is  a  complex  signal.    This  is  because 

•  The  WVD  satisfies  the  following  symmetry  relations: 

W(t,u)  =  W(t,-u)    if    s(t)     is  real,  i.e.,  if    5(w)  =  S(-w), 
and 

W{t,w)  =  W(-t,u)     if    S(u)     is  real,  i.e.,  if     s(t)  =  s(-t). 

•  The  marginal  property  is  satisfied  in  both  time  and  frequency.  Hence  the  total 
energy  property  is  also  satisfied. 

•  The  time  and  frequency  shift  invariance  properties  are  satisfied. 
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For  s(t)  =  A(t)e?*w  the  conditional  spread  in  frequency  is  given  by 

-\2 


°l,t  =  ("2)t  -  l(u)tY  =  0.5 


A'{t) 

A(t) 


A"(t) 

A(t) 


(11.30) 


where  "/"  denotes  the  time  derivative.  Since  Equation  11-30  may  result  in  neg- 
ative values,  it  cannot  be  properly  used  as  a  measure  of  bandwidth.  Therefore, 
while  WVD  can  be  used  to  compute  a  reliable  mean  frequency,  it  cannot  be 
used  to  obtain  reliable  values  for  the  spread  of  those  frequencies. 

•  WVD  of  sum  of  two  signals:  If  s(t)  =  Si(t)  +  s2{t),  then  the  WVD  is  given  by 
W(t,u>)  =  Wu(t,u)  +  W»(*,w)  +  W12(t,u)  +  W2l(t,u)),  (11.31) 

where 

t-T\  (t  +  T* 


w^)=\h\iy 


,-JWT 


dr 


(11.32) 


is  the  cross  WVD.  The  cross  WVD  is  complex- valued,  but  Wu  =  W2\.  Hence, 
W\2{t,w)  +  W2\{t,u))  is  real  and  we  can  express  W(t,u)  as: 


W(t,u)  =  Wu(t,u)  +  Wn&u)  +  2Re  [Wl2{t,u)] . 


(11.33) 


The  additional  term  2Re  [^2(^,0;)]  is  often  called  the  cross  term  or  the  inter- 
ference term. 

Note:  The  signal  is  said  to  be  composed  of  more  than  one  signal  (i.e.,  multi- 
component  signal)  if  it  can  be  segmented  into  separate  regions  in  the  time-frequency 
plane. 

Note:  The  WVD  at  time  t  depends  on  the  values  of  the  signal  at  the  time 
moments  {t  —  r/2)  and  (t  +  r/2).  Even  if  there  is  no  noise  at  time  t,  the  WVD  will 
reflect  the  noise  component  existing  at  those  time  moments  (t  —  r/2)  and  (t  +  r/2). 
There  are  other  types  of  TFDs,  all  forming  one  class,  called  the  general  Cohen's 
class.  Cohen's  class  incorporates  different  kinds  of  TFDs  which  differ  only  in  the 
definition  of  the  window  function.  Window  functions  may  be  chosen  to  optimize  the 
TFD  according  to  certain  criteria.  One  example  of  this  criteria  might  be  to  minimize 
cross  terms. 
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C.      WAVELET  ANALYSIS 

Wavelet  analysis  is  a  new  trend  for  representing  non-stationary  signals.  It  is 
widely  applied  for  many  areas  of  signal  processing.  In  the  following  discussion  we 
introduce  the  basic  concepts  of  wavelet  analysis. 

1.       Frames 

A  basis  for  a  vector  space  V  is  a  set  S  of  vectors  in  V  such  that  S  is  linearly 
independent  and  S  spans  V  [Ref.  41].  Therefore,  a  basis  can  express  any  vector 
in  the  space  V  by  a  linear  expansion.  The  resultant  set  of  expansion  coefficients  is 
unique.  Completeness  and  uniqueness  of  the  expansion  coefficients  are  directly  related 
to  the  exact  reconstruction  of  the  original  vector  from  the  expansion  coefficients. 
This  property  is  required  in  signal  processing  applications  which  use  the  expansion 
algorithm  to  reconstruct  the  original  signals.  For  simple  computation  of  the  expansion 
coefficients,  the  vectors  of  S  are  required  to  be  orthogonal.  Orthogonality  will  make  it 
possible  to  compute  the  expansion  coefficients  (analysis)  using  simple  inner  products 
with  the  basis  vectors.  Reconstruction  (synthesis)  and  decomposition  (analysis)  are 
done  using  the  same  set  of  basis  vectors. 

Orthogonality  of  the  basis  vectors  (or  basis  functions)  is  defined  as  follows: 
the  set  of  basis  vectors  Vk(t)  is  said  to  be  orthogonal  if 

{vk(t),  vm(t))  =  ak,m6(k  -  m),  (11.34) 

where  k,  m  are  integers,  ak,m  is  a  constant  and  S  is  the  Kronecker  delta  function.  The 
delta  function  6(k  —  m)  will  be  zero  everywhere  except  where  m  =  k  it  will  be  one. 
Since  the  orthogonality  is  not  a  condition  for  having  a  set  of  basis  vectors,  we  may 
find  a  set  of  non-orthogonal  basis  vectors.  For  any  orthogonal  set  of  basis  vectors  Vk(t) 
the  finite  energy  time  signal  f(t)  can  be  expanded  using  the  set  of  basis  functions 
vk(t)  as: 
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where 

o*-</(*W*)>, 

and  ||vjfc(£)||  is  the  1?  norm  of  the  vector  vk(t).  Further,  if  we  assume  unit  L2  norm  for 
the  vectors  Vk(t),  the  vectors  become  an  orthonormal  set,  and  the  expansion  formula 
will  be 

/(*)  =  £  <*«*(*),  (H-36) 

k 

where 

ak  =  (f(t),vk{t)). 

It  is  worth  mentioning  here  that  due  to  orthogonality  (or  orthonormality)  both 
the  analysis  and  synthesis  of  the  expansion  are  carried  out  using  the  same  set  of  basis 
vectors. 

For  a  biorthogonal  basis,  we  need  two  different  sets  of  basis  vectors,  the  original 
set  Vk(t)  (the  basis)  for  the  expansion  and  a  dual  set  Vk(t)  (dual  basis)  for  reconstruc- 
tion (or  synthesis).  The  sets  of  the  basis  and  the  dual  basis  are  not  orthogonal,  but 
each  vector  and  its  dual  are  orthogonal.  The  expansion  is  given  by 

/(t)  =  Efl*v»W.  (IL37) 

k 

where 

ak  =  (f{t),vk{t)). 

If  the  set  S  of  the  expansion  vectors  does  not  satisfy  the  linear  independence 
condition  stated  above,  then  the  set  is  called  a  "frame."  Expansion  with  frames  does 
not  maintain  uniqueness  of  expansion  coefficients  [Ref.  37].  Also,  the  exactness  of 
the  reconstruction  will  be  replaced  with  an  approximate  representation.  For  more 
details  about  the  theory  of  frames  see  [Ref.  42,  43].  Frames  do  not  satisfy  the 
Parseval's  theorem  of  energy,  which  states  that  the  energy  in  the  original  function 
will  be  distributed  among  the  expansion  coefficients  without  loss  or  gain.  For  frames, 
the  total  energy  in  the  expansion  coefficients  will  have  two  bounds  A  and  B  such 
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that: 

^ll/ll2  <£!</(*),  M*))l2<£||/lf  (H.38) 

k 

where  0  <  A,  B  <  oo.  Note  that  if  A  =  B,  then  the  frame  is  called  a  tight  frame  and 

^ll/ll2  =  £K/(*W')>l2-  (H.39) 

k 

If  A  =  B  =  1,  then  the  frame  becomes  an  orthonormal  set  and 

ll/ll2  =  £l(/(*W*)>l2,  (H.40) 

k 

which  is  the  the  Parseval's  theorem. 

There  are  different  categories  of  wavelets,  orthogonal  wavelets,  non-orthogonal, 
and  biorthogonal  wavelets.  The  Daubechies  family,  Symmlet,  Coiflet,  and  Meyer 
wavelets  are  examples  of  orthogonal  wavelets,  while  the  Morlet  wavelet  is  an  example 
of  a  non-orthogonal  wavelet.  As  mentioned  earlier,  not  all  applications  of  signal 
expansions  require  perfect  reconstruction.  Applications  such  as  signal  identification 
and  classification  are  examples  of  this  type.  Therefore,  unless  the  orthogonality  is 
required  to  meet  another  criteria,  we  can  exploit  non-orthogonal  wavelet  transforms. 

2.       Continuous  Wavelet  Transform 

The  continuous  wavelet  transform  (CWT)  forms  the  mathematical  basis  for 
wavelet  analysis.  The  concept  behind  wavelet  analysis  is  that  all  basis  functions  can 
be  generated  from  a  single  function  called  the  mother  wavelet,  usually  denoted  by 
ip(t).  Other  wavelets  can  be  generated  using  two  distinct  operations;  scaling  and 
translation.  Scaling  means  expanding  (dilating)  or  compressing  the  wavelet  function 
according  to  a  specific  scaling  value.  The  scale  will  be  denoted  by  s.  The  translation 
allows  shifting  of  the  (scaled)  wavelet  to  a  desired  position  in  time.  This  shift  will  be 
denoted  by  k.  The  scaled  and  translated  wavelet  is  denoted  by 

«*»(*)  =  -^  fcjty  ,  (11.41) 

where  >/s  is  a  normalization  factor  to  maintain  the  L2  norm  of  ips,k{t)  to  be  constant 
at  any  scale  s  . 
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The  integral  form  of  CWT  of  the  finite  energy  signal  f(t)  with  respect  to  the 
wavelet  function  if)(t)  is  given  by  [Ref.  42,  46] 

/oo 
f(t)ips,k(t)dt,  (11.42) 

-00 

or 

1         /-OO  fi  _  h\ 

Wf(s,  k)  =  ■—  ^  /(t)tf  (— ^J  eft.  (11.43) 

The  wavelet  analysis  is  carried  out  by  computing  inner  products  (projections) 
between  the  signal  and  the  wavelet  functions.  We  can  also  interpret  the  wavelet  anal- 
ysis as  a  linear  operation  which  transforms  the  signal  using  modified  kernel  functions, 
where  the  kernel  of  the  transform  is  the  mother  wavelet,  and  the  modifications  are  the 
scaling  and  the  translation  operation  [Ref.  44].  If  the  wavelet  function  satisfies  the 
admissibility  condition,  we  can  apply  the  rule  of  the  resolution  of  identity  to  recover 
the  signal  from  its  wavelet  transform.  The  admissibility  condition  is  given  by  [Ref. 

44] 

r°°  \ty(u)\2 
C*  =  /      '    ,    ,     dw  <  oo,  (11.44) 

./-oo        M 


where  &(w)  is  the  Fourier  transform  of  ip(t).  This  condition  will  be  true  if  ^(0)  =  0 
or  equivalently 

/OO 
ip(t)dt  =  0.  (11.45) 

-oo 

The  admissibility  condition  implies  that  the  wavelet  has  a  zero  mean  (i.e.,  DC  com- 
ponent is  zero).  Therefore,  it  can  be  interpreted  as  a  bandpass  function.  This  implies 
that  the  wavelet  must  be  an  oscillatory  function.  The  recovery  (synthesis)  formula, 
extracted  from  the  resolution  of  identity  rule,  is  given  by 

/•OO      /"OO 

/(*)=/      /     Wf(s,k)ips,k(t)dsdk.  (11.46) 

•Zoo    -Zoo 

3.       Properties  of  Continuous  Wavelet  Transform 

The  following  properties  of  the  wavelet  ransform  are  given  below  without  proof. 
More  details  can  be  found  in  [Ref.  45]. 
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1.  Linearity:  The  wavelet  transform  is  a  linear  operation  on  the  analyzed  signal. 
If  the  signal  f(t)  is  given  as: 

/(t)  =  /i(t) +  /»(*),  (IL47) 

then  the  wavelet  transform  of  f(t)  is  given  by 

Wf(s,  k)  =  Wfl(s,  k)  +  Wh(s,  k).  (11.48) 

2.  Shift-invariance:  The  CWT  is  shift  invariant.  If  g(t)  =  f(t  —  t0),  then  it  has 
a  CWT  given  by 

Wg(s,k)=Wf(s,k-t0).  (11.49) 

3.  Time-scaling-invariance:  If  g(t)  =  ^/(~)3  then  it  has  a  CWT  given  by 

Wg(s,k)  =  Wf  (-,-).  (11.50) 


4.  Parseval's  relation  for  the  CWT  becomes: 

i    /-00  r30 

implying  that  the  CWT  conserves  the  signal  energy. 


roo  1         roo     roo  riddle 

/    \\f(t)\\2dt  =  ^r        /    \\Wf(s,k)\\2^f-  (H-51) 


5.  Time  localization:  From  the  definition  of  CWT  and  by  the  sifting  property  of 
the  Dirac  delta  function,  the  CWT  of  an  impulse  at  time  t0  will  be  [Ref.  47] 

i  ^klj.  This  means  that  the  response  of  the  CWT,  at  any  scale  s,  will  be 
a  scaled  and  time  reversed  replica  of  the  mother  wavelet  centered  at  location 
t0  on  the  shift  axis  (k).  Therefore,  the  ability  to  define  the  location  t0  will 
depend  on  the  width  of  the  scale  wavelet.  And  since,  a  smaller  scale  value  will 
result  in  a  shorter  wavelet,  the  time  localization  will  be  more  precise. 

6.  Frequency  localization:  In  fact,  the  uncertainty  principle  controls  the  time 
and  frequency  resolution  of  the  wavelet  transform  at  any  scale,  as  the  smaller 
the  scale  is  the  shorter  the  wavelet  and  the  wider  its  frequency  representation 
are.  Consequently,  the  shorter  wavelet  has  worse  frequency  localization  than 
the  original  wavelet.  Generally,  the  time  localization  is  better  at  smaller  scale 
values  while  the  frequency  localization  is  better  at  larger  scale  values. 

4.       Scalogram 

Using  a  scalogram,  the  finite  energy  signal  f(t)  is  represented  by  the  distribu- 
tion of  |W/(s,  A;)|2  over  the  time-scale  plane.  From  Parseval's  theorem  for  the  wavelet 
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transform,  Equation  11-51,  integration  of  |W)(s,  k)\2  with  respect  to  the  differential 
^r  provides  the  total  energy.  By  the  definition  of  the  scaled  wavelet  function  in 
Equation  11.41,  the  scale  value  s  is  inversely  proportional  to  the  frequency.  Hence, 
the  differential  ^f  is  proportional  to  the  differential  of  the  frequency.  Therefore,  the 
quantity  |W/(s,  k)\2  may  be  viewed  as  a  spectral  density  in  units  of  power  per  Hertz 
[Ref.  44].  Consequently,  the  scalogram  represents  the  power  spectral  density  of  the 
signal  over  the  time-scale  plane.  Finally,  we  have  the  quantity 

1        /-oo  /7c 

rr       IIW7(».*)IIJ3  (IL52) 

G,/,  Joo  Sz 

which  represents  the  signal  instantaneous  power,  while  the  quantity 

1  roo     roo 

rTS        /     \\Wf(s,k)\\2dk  (11.53) 

O^S     Joo    Joo 

represents  the  portion  of  the  signal  energy  contained  within  the  scale  s.  This  fact 
is  exploited  in  identifying  the  scale  of  each  hop.  In  summary,  the  CWT  is  a  linear, 
time-shift-invariant,  time-scaling-invariant,  and  frequency-scaling-invariant  operator. 

5.       Discrete  Wavelet  Transform 

The  CWT  is  defined  by  an  integral  transform  over  continuous  variables  in 
the  scale  s  and  the  time  shift  k.  Hence,  performing  the  CWT  requires  expensive 
computations.  Therefore,  in  practice  a  discrete  grid  for  s  and  k  is  used.  A  widely 
accepted  discretization  is  to  specify  s  =  s™  and  k  =  nk0s™,  where  m  and  n  are 
integers,  and  s0  >  1,  and  k0  >  0  [Ref.  45].  Furthermore,  if  we  select  s0  =  2  and 
k0  =  1  we  obtain  the  well-known  dyadic  wavelet  sampling  (tiling)  grid.  Hence  the 
scaled  and  translated  wavelet  indexed  by  m  and  n  is  given  by 

V>TO,„(t)  =  2=?iP{2-mt  -  n).  (11.54) 

Figure  1  shows  the  time-frequency  tiling  of  the  STFT  while  Figure  2  shows  the  time- 
scale  tiling  of  the  wavelet  transform. 
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In  digital  signal  processing  there  are  two  related  concepts  that  help  under- 
standing the  wavelet  transforms:  the  multiresolution  analysis  and  the  filter-bank 
theory.  These  concepts  are  discussed  next. 

6.       Multi-Resolution  Concepts  in  Discrete  Wavelet 
Transforms 

Introduction.  We  will  briefly  introduce  the  WT  from  the  perspective  of  the  mul- 
tiresolution analysis  (MR).  The  signal  (or  the  time  function)  f(t)  is  expanded  in 
terms  of  the  wavelet  functions.  These  wavelets  have  a  frequency  bandpass  shape,  so 
they  result  in  a  set  of  successive  details  of  the  signal.  For  the  approximation  we  need 
a  special  basis,  called  the  scaling  function  <j>(t),  which  is  not  a  wavelet.  It  has  a  low 
pass  frequency  behavior  and  performs  averaging.  The  discrete  wavelet  transform  is 
given  by 

OO  OO  00 

/(*)=    £   c(k)fa(t)  +  Y,   £   d(j,k)i>{j,k)(t),  (11.55) 

k=—oo  j=0  k=—oo 

where  j  and  k  are  integers,  the  coefficient  c(k)  constitutes  the  coefficients  of  the 
approximation,  while  d(j,  k)  constitutes  the  coefficients  of  the  added  details  or  equiv- 
alent^ the  fine  resolutions  [Ref.  47].  If  the  wavelets  and  the  associated  scaling 
functions  form  an  orthonormal  set  of  basis  functions  the  coefficients  are  given  by 

c(fc)  =  (/(*).  &(*)>.  (11.56) 

and 

d(j,k)  =  {f(t),1>jtk(t)).  (11.57) 

The  expansion  form  of  the  DWT  is  related  to  the  subspace  representation.  Let 
us  define  the  following  subspaces:  Vj  is  the  scaling  space  at  the  jth  level,  and  Wj  is  the 
wavelet  space  at  thejth  level.  Let  Vj+i  =  Vj+Wj.  This  means  that  the  approximation 
at  the  (j  +  l)st  level  (or  scale)  can  be  represented  by  a  coarser  approximation  Vj  and 
a  coarser  detail  Wj,  both  at  the  jth  level.  Note  that  the  subspaces  Vj  are  spanned  by 
the  scaling  functions  fa  (t)  and  the  subspaces  Wj  are  spanned  by  the  wavelet  functions 

iM*). 

22 


Subspaces  have  to  satisfy  a  set  of  conditions  to  become  an  MR  representation. 
These  conditions  are  [Ref.  46]: 

1.  Orthogonality  and  completeness: 

VjHWj    =    {0}. 
VjQWj   =   vj+1. 

00 

V}®Y^Wj    =    L2.  (11.58) 

j=0 

2.  Scale-invariance:  If  /(*)  €  Vj  then  f(2t)  £  Vj+l. 

3.  Shift-invariance:  If  f(t)  €  Vj  then  f(t  —  k0)  €  Vj  for  a  given  constant  k0  . 

Seeding  and  Wavelet  Equations.  The  MR  subspace  representation  leads  to  a 
method  to  formulate  two  equations  in  terms  of  the  unknown  scaling  and  wavelet 
functions.  From  the  orthogonality  and  completeness  property  we  deduce  that  Vo  C 
V\.  This  means  any  function  living  in  V0,  e.g.,  (j>o(t),  can  be  expressed  as  a  linear 
combination  of  functions  living  in  Vi,  e.g.,  <f>\(t).  Hence,  by  substitution  we  get 

(f>(t)  =  V2J2  h0(h)(f>{2t  -  fc),  (11.59) 

k 

where  (f>o(t)  =  <f>(t)  and  <f>\{t)  =  <f>{2t).  This  equation  includes  two  different  scales  of 

the  scaling  function  and  is  known  by  many  names:  the  scaling  equation,  the  dilation 

equation,  the  refinement  equation,  or  the  multiresolution-analysis  equation.     The 

coefficients  h0(n)  are  called  the  scaling  filter  coefficients.  Using  the  condition,  V0  © 

Wo  =  Vi,  which  implies  that  Wq  C  Vi,  we  have 

if>(t)  =y/2j2  h^k^t  -  k)  (11.60) 

k 

also  is  called  the  wavelet  equation.  The  coefficients  h\(n)  are  called  the  wavelet  filter 
coefficients. 

In  summary,  an  MR  approach  leads  to  a  complete  orthonormal  set  of  the 
scaling  and  wavelet  functions.  Moreover,  if  the  scaling  function  has  compact  support, 
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the  wavelet  function  will  be  built  up  of  a  finite  number  of  the  scaling  multipliers. 
Therefore,  the  scaling  and  the  wavelet  functions  can  be  realized  by  finite-impulse- 
response  (FIR)  filters.  For  more  details  about  wavelet  and  scaling  equations  see  [Ref. 
47,  42]. 

For  discrete  data  the  filter-bank  concept  leads  to  a  simple  method  for  com- 
puting the  wavelet  coefficients.  The  wavelet  function  is  replaced  by  the  coefficients 
of  the  wavelet  filter  hi(n)  and  the  scaling  function  is  replaced  by  the  coefficients  of 
the  scaling  filter  h0(n).  In  [Ref.  46,  47,  42],  the  following  two  recursive  equations  are 
obtained: 

c(j> k)    =    ^2h0(m-  2k)c(j  +  1,  m) 

m 

d(j,k)    =    £Mm-2fe)c(j  +  l,m).  (H.61) 

m 

These  two  recursive  equations  enable  us  to  compute  the  jth  scale  wavelet  trans- 
form from  the  (j  +  l)st  scale  using  ho(m)  and  hi(m).  The  factor  2  in  front  of  the 
parameter  k  decimates  the  output  wavelet  coefficients  by  2.  Figure  3  shows  the  re- 
alization of  the  recursive  equations.  This  realization  is  equivalent  to  the  well-known 
2-band  analysis  bank  of  filter  implementation.  The  implementation  of  the  discrete 
wavelet  transform  using  the  recursive  procedure  is  known  as  Mallat's  algorithm.  Fig- 
ure 4  shows  a  multi-stage  wavelet  analysis  bank  and  the  corresponding  ideal  subspace 
designation  (or  spectral  decomposition). 

7.       Daubechies  Wavelet  Family 

The  Daubechies  wavelet  family  is  a  compactly  supported  orthonormal  set  of 
wavelet  functions  [Ref.  42].  The  Daubechies  wavelet  are  generated  by  solving  the 
scaling  and  wavelet  equations.  An  additional  set  of  constraints  is  applied  to  satisfy  the 
maximum  number  of  vanishing  moments  for  each  wavelet.  Each  Daubechies  wavelet 
is  assigned  a  number  related  to  the  number  of  vanishing  moments,  hence,  each  wavelet 
is  denoted  as  Daub-N.  There  are  two  conventions  in  assigning  the  number  N;  the  first 
one  is  to  let  iV  refer  to  the  number  of  vanishing  moments;  the  second  one  is  to  let  N 
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refer  to  the  length  of  the  wavelet  filter  (i.e.,  the  number  of  wavelet  coefficients).  We 
will  adopt  the  first  designation  since  it  is  adopted  by  the  MATLAB  Wavelet  Toolbox, 
where  each  Daubechies  wavelet  is  expressed  by  the  coefficients  of  its  scaling  filter. 
The  coefficients  of  the  wavelet  filters  are  the  interpolated  high  pass  versions  of  the 
scaling  filter.  The  Daubechies  wavelet  of  order  N  has  2N  coefficients,  and  has  finite 
support  over  [0,  27V  —  1].  The  number  of  vanishing  moments  is  an  indication  of  the 
smoothness  of  the  the  wavelet  filter,  since  the  number  of  vanishing  moments  will 
imply  the  number  of  zeros  of  ^(lj)  at  u  =  tt.  The  higher  the  order  the  longer  and 
smoother  is  the  Daubechies  filter. 
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Figure  1.  Time-Frequency  tiling  of  the  STFT. 
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Figure  2.  Time-Scale  tiling  of  the  wavelet  transform. 
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Figure  3.  Recursive  computation  of  wavelet  coefficients. 
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Figure  4.  Multiple-Stage  wavelet  analysis  bank  and  its  subspace  designation. 
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III.         WAVELET  TRANSFORMS  AND 
CORRELATION  FUNCTIONS 

Correlating  two  functions  results  in  a  measure  of  similarity  between  them 
and  is  accomplished  by  evaluating  an  inner  product.  Correlation  analysis  has  been 
used  extensively  in  the  signal  processing  and  communication  area,  e.g.,  in  spectral 
estimation,  system  modeling,  signal  detection  and  parameter  estimation. 

The  Wiener-Khinchin  theorem  relates  the  signal's  auto-correlation  function 
and  power  spectral  density  for  a  stationary  stochastic  process,  through  a  Fourier 
transform.  Due  to  their  time-varying  spectra,  FH  signals  are  non-stationary  and  can 
be  represented  using  time-frequency  distributions.  Traditionally,  Fourier  kernels  are 
used  in  the  time-frequency  distributions. 

Wavelet  decomposition  can  be  used  to  represent  non-stationary  signals  over  the 
time-scale  plane.  Therefore,  we  will  consider  the  wavelet  transform  of  the  correlation 
function  as  an  alternative  for  non-stationary  signal  representation. 

In  this  chapter,  we  will  introduce  different  choices  for  the  correlation  function 
and  derive  the  corresponding  wavelet  transform  response.  In  the  derivation  of  the 
wavelet  transform  we  will  use  the  analytic  definition  of  the  correlation  function  and 
the  wavelet  function  without  specifying  a  functional  form. 

A.      CORRELATION  FUNCTIONS 

1.       Definitions  of  Auto-Correlation  Functions 

Depending  on  the  underlying  process,  various  definitions  are  given  to  the 
auto-correlation  function  (ACF).  The  process  may  be  deterministic,  stochastic,  finite- 
energy,  infinite-energy,  non-time- varying  (stationary)  or  time- varying  (non-stationary). 
Next,  we  will  introduce  the  different  definitions  of  the  ACF.  The  main  reference  for 
this  topic  is  [Ref.  51]. 

1.  Deterministic  ACF:  The  ACF  is  defined  for  deterministic,  infinite-energy,  and 
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non-time- varying  signals  by 

R{T)  =  £%o  ?T  LTX(t)7r{t  +  T)dt  '  (HL1) 

where  "T"  is  the  observation  interval  and  "r"  is  the  time  lag  or  the  time  delay. 
For  finite  energy  signals  the  ACF  is  defined  as: 

x(t)x*{t  +  r)dt.  {111.2) 

-00 

For  stationary  signals  or  processes  (i.e.,  wide-sense  stationary)  R(t,r)  is  a 
function  of  the  time  delay  only,  hence,  it  is  denoted  by  R(t).  For  a  non- 
stationary  signal,  R(t,  r)  will  be  represented  by  a  2-D  surface  over  the  time 
"t"  and  the  time  delay  "r". 

2.  Stochastic  ACF:  The  ACF  of  a  stochastic  process  is  the  correlation  between 
two  samples  of  the  process  taken  at  t\  and  t2,  and  is  defined  as: 

R{t1,t2)  =  E{x(t1)x*(t2)},  (111.3) 

where  E  is  the  expectation  operator  and  "*"  stands  for  the  complex  conju- 
gation. For  a  stationary  (or  wide-sense  stationary)  process,  R(ti,t2)  depends 
only  on  the  time  lag  r  =  t2  —  ti,  resulting  in  a  stationary  spectrum.  For  a  non- 
stationary  process,  R{ti,t2)  depends  on  the  time  instants  t\  and  t2,  resulting 
in  a  time- varying  spectrum. 

2.       Properties  of  Auto- Correlation  Functions 

We  will  briefly  introduce  the  significant  properties  of  the  ACF  for  finite-energy 
signals.  For  more  details  see  [Ref.  36,  52]. 

1.  Conjugate  symmetry: 

R(t,r)  =  R{t,-T)*,  .    (III.4) 

i.e.,  the  real  part  of  R(t,  r)  is  an  even  function  with  respect  to  the  lag  r  while 
the  imaginary  part  is  an  odd  function. 

2.  Mean-squared  value: 

R(t,r)\T=0  =  E{\x(t)\2}.  (III.5) 

3.  Positive  definiteness:   The  ACF,  denoted  by  R(ti,t2),  is  said  to  be  positive 
definite  if 

Y^aia^Riti.t^yO  (III.6) 

for  any  ai,a,j  ^  0.  Positive  definite  property  of  the  ACF  guarantees  that  the 
spectrum  of  the  signal  is  non-negative. 
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4.  Wiener-Khinchin  theorem: 

/oo 
R(r)e-j2^dr.  (III.7) 

-00 

The  power  spectral  density  function  Sxx(f)  of  the  stationary  signal  is  ob- 
tained by  Fourier  transforming  its  (positive  definite)  ACF.  Note  that  the  non- 
stationary  signals  have  time- varying  spectra,  therefore,  no  power  spectral  den- 
sity is  defined. 

3.       The  Instantaneous  Correlation  Function 

The  ACF  of  deterministic  and  stochastic  processes  are  computed  using  time 
domain  averaging  and  the  expectation  operator,  respectively.  This  means  that  a 
smoothing  process  has  to  be  applied  to  compute  the  correlation  functions. 

The  instantaneous  correlation  function  (ICF)  does  not  use  an  averaging  oper- 
ation (integration  nor  expectation).  The  instantaneous  correlation  function  is  defined 
as  the  product  of  two  samples  of  the  process.  These  two  samples  are  drawn  at  two 
time  instants  centered  about  time  t.  The  instantaneous  correlation  function  R%  (t,  r) 
is  defined  as: 

#feT)  =s(t  +  j)  *•(*-£),  (IH.8) 

where  "i"  stands  for  the  instantaneous  nature  of  the  correlation  function. 
Rl{t,r)  satisfies  the  following  properties: 

1.  Conjugate  symmetry,  i.e., 

Ri(t,T)  =  Ri*(t,-r).  (III.9) 

2.  Squared  value, 

Ri(t,r)\T=:0  =  \x(t)\2,  (111.10) 

i.e.,  at  r  =  0  ,  Rl(t,r)  will  represent  the  instantaneous  power  of  the  signal  at 


«j.» 


time  "t 

If  x(t)  is  a  sinusoidal  signal  then  multiplication  of  the  instantaneous  values  of 

x(t)  in  the  sense  of  Equation  III-8  will  generate  cross  terms  in  the  ICF.  For  example, 

the  real-valued  sinusoidal  signal  x(t)  =  Acos(ut)  has  an  ACF  given  by: 

A2 
R(t)  =  —cos{ut),  (III.ll) 
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while  the  ICF  is  given  by: 

A2 
R*(t,  r)  =  —  [cos(2ut)  +  cos(o;r)] .  (111.12) 

The  ACF  of  a  single  sinusoidal  signal  has  only  one  component  and  no  cross 
term,  while  the  ICF  has  cross  terms.  If  the  signal  x(t)  is  represented  by  its  analytic 
form,  say  x{t)  =  Aexpjcut,  then  its  ICF  is  given  by: 

Rl(t,  r)  =  A2  exp  (jut).  (111.13) 

That  is,  the  ICF  of  a  single  complex  exponential  signal  has  no  cross  term.  Therefore, 
to  minimize  cross  terms  it  is  recommended  to  use  the  analytic  form  of  the  input 
process  in  the  computation  of  the  Time- Frequency  Distributions  [Ref.  49].  In  this 
dissertation,  we  will  focus  on  the  analytic  form  of  the  signals.  Therefore,  the  first 
step  of  the  processing  scheme  will  be  transforming  the  intercepted  (real)  signal  using 
the  Hilbert  transform  technique. 

Another  point  worth  mentioning  is  the  exploitation  of  the  conjugate-symmetry 
property  of  the  ICF.  In  computing  the  surface  of  the  ICF  we  need  to  compute  only 
half  the  ICF  surface,  corresponding  to  positive  (or  negative)  time  lag  V.  The  other 
half  can  be  generated  by  the  conjugate-symmetry  property.  This  implies  that  the 
second  half  of  the  surface  of  the  ICF  carries  no  additional  information.  Therefore, 
the  processing  scheme  may  be  applied  to  only  one  half  of  the  ICF  surface  which  will 
save  computational  cost  and  storage. 

B.      WAVELET  TRANSFORMS  OF  CORRELATION 
FUNCTIONS 

In  Chapter  II,  we  discussed  the  linearity  property  of  the  wavelet  transform 

which  allows  the  application  of  linear  system  theory  to  the  wavelet  transform.  In  this 

section,  we  will  investigate  the  wavelet  transform  for  the  different  definitions  of  the 

correlation  function. 
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1.        The   Wavelet    Transform   of  the   Auto-Correlation 
Function 

The  wavelet  transform  of  the  stochastic  ACF  of  a  stationary  process  is  ad- 
dressed in  [Ref.  48],  similarly,  we  will  address  the  ACF  for  a  deterministic  signal. 
Wavelet  transform  of  the  ACF  will  have  similar  results.  For  a  deterministic,  finite- 
energy,  stationary  signal,  the  (time-averaged)  ACF  is  given  by: 

/oo 
x{t)x*{t  +  r)dt.  (111.14) 

-oo 

From  the  Wiener-Khinchin  theorem  we  get 

/oo 
R(r)e-j2*fTdT.  (III.  15) 

-oo 

Let  Wxx(s,£)  denote  the  Wavelet  transform  of  R(j).  Note,  the  subscript  "xx" 
in  Wxx(s, £)  stands  for  the  wavelet  transform  of  the  ACF  of  x(t)  (in  contrast  to 
Wx(s,£)  which  denotes  the  wavelet  transform  of  x(t)).  The  wavelet  basis  function 
is  denoted  by  g(r).  The  wavelet  transformation  will  be  carried  over  the  time  lag 
variable  "r"  to  the  wavelet  shift  variable  "f  and  the  wavelet  scale  "s".  Wxx(s,£)  is 
then  given  by: 

Wxx(s,  £)  =  -~  jT  R(T)g*  ( Lzij  dr .  (111.16) 

Let  G(f)  denote  the  FT  of  g(r),  then 

FT{g  (L^)>  =  1*1  G(sf)e-^*.  (111.17) 

For  positive  scale  values  ,"s",  g  f2^)  is  given  by: 

g  (j^J  =  f°  sG{sf)e-j2*feej2*fTdf .  (III.  18) 

We  can  write  Wxx(s,  £)  as: 

/OO        /-OO 
/      R(T)e-j2*fTdrG*  {sf)ej2"fTdf,  (III.  19) 

-00  J— 00 


or 


/oo 
sxxU)G\sfy2*ftdf .  (111.20) 

-00 


35 


Equation  III. 20  has  the  form  of  an  inverse  FT  from  the  variable  /  back  to  the 
variable  t.  Therefore,  we  can  write 

Wxx(s,£)  =  FT-'i^sS^G^sf))}  .  (111.21) 

From  Equation  III. 21  we  deduce  the  following: 

1.  The  wavelet  transform,  of  a  finite-energy  signal,  at  any  wavelet  scale  s  ^  0, 
represents  a  linear  filtering  operation  using  a  band  pass  filter  whose  impulse 
response  is  the  (time-reversed)  wavelet  function  at  the  scale  s.  Equivalently, 
the  filter  has  a  frequency  response  given  by  the  FT  of  the  scaled  wavelet. 

2.  The  wavelet  transform  of  the  ACF,  R(r),  of  the  stationary  finite-energy  signal 
x(t),  gives  a  band  pass  filtered  version  of  the  power  spectral  density  Sxx(f)  of 
this  signal  (up  to  a  constant,  y/s),  the  used  band  pass  filter  is  dependent  on 
both  the  chosen  wavelet  function  and  the  chosen  wavelet  scale. 

2.       The  Wavelet  Transform  of  the  Instantaneous 
Correlation  Function 

The  Wigner-Ville  Distribution  (WVD)  is  used  to  represent  non-stationary  pro- 
cesses since  they  are  characterized  by  their  time- varying  spectra.  The  WVD  applies  a 
one-dimensional  Fourier  transformation  to  the  ICF.  The  Fourier  transform  is  carried 
out  taking  the  time  delay  r  to  the  frequency  /,  keeping  the  global  time  t  unchanged. 
This  allows  for  displaying  the  time  evolution  of  the  spectrum  of-  the  signal.  For 
one-dimensional  time  signals,  the  one-dimensional  wavelet  transform  carries  out  a 
transformation  from  one  global  time  variable  Hn  to  the  two  wavelet  variables,  the 
shift  i  and  the  scale  s.  Consequently,  the  signal  is  represented  by  a  time-scale  distri- 
bution in  the  wavelet  domain.  The  wavelet  domain  is  called  the  time-scale  domain. 
For  the  two-dimensional  surface,  indexed  by  time  "£"  and  the  delay  "r" ,  we  carry  out 
the  wavelet  transformation  along  the  time  delay.  This  permits  a  display  as  a  func- 
tion of  time.  In  this  section,  we  will  derive  the  relation  between  the  time-frequency 
representation  of  the  signal  x{t)  using  the  WVD  and  the  time-scale  representation  of 
the  same  signal  using  the  wavelet  transform  of  the  ICF. 
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Let  Vx(t,  f)  denote  the  WVD  of  the  signal  x(t), 

Vx{t>  f)  =  Co  *((*  "  l)  *'  (*  +  5)  e"J27r/TdT-  (nL22) 

Since  the  ICF  of  the  signal  x(t)  is  defined  as: 

(111.23) 
where  r  is  time  shift  (or  time  delay)  relative  to  time  t.  The  WVD  may  be  defined  as: 

/oo 
Rx{t,r)e-^dT.  (111.24) 

-00 

The  wavelet  function  g{r),  has  been  scaled  by  s  and  shifted  by  I  along  the  r  axis,  to 
form  gs,e{T)  which  is  given  by: 

9sAt)  =  Tsg  fr^ ' 

Let  G(f)  denote  the  FT  of  g(r),  thus 


for  s  >  0  . 


(111.25) 


Then 


FT{g  l^)}  =  s  G(sf)e-^e 


L-l\  =  j°°  sG(sf)e-j2*feej2*fTdf 


(111.26) 


(111.27) 


Let  Wxx(t]s,£)  be  the  wavelet  transform  of  the  ICF  Rx(t,  r).   Note  that  the 
superscript  "z"  indicates  the  instantaneous  feature  of  the  ICF.  W*x(t\  s,  £)  is  given  by: 


WLfrsJ)  =  jlf"  Rx{t,T)g*  f^-J  dr 


(111.28) 


yj S  J—  oo  V        S 

Substituting  p*  f2"^)  fr°m  Equation  III. 27  into  Equation  III. 28  and  exchang- 
ing the  integration  operations  we  get 

/oo  r  /•oo  "I 

G*(sf)    /      Rx(t,r)e-j2"fTdT  ^f£4f.  (111.29) 

-00  ly— 00 

Substituting  Equation  111.24  into  Equation  III. 29  we  have 

/oo 
G*(sf)Vx(t,f)e^2^edf.  (111.30) 

-00 
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Equation  III. 30  is  in  the  form  of  an  inverse  Fourier  transform.  Therefore,  W*x(t]  s,£) 
and  y/sG*(sf)Vx(t,f)  are  a  Fourier  transform  pair  with  respect  to  the  variable  u£" , 
i.e.,  af.  This  relation  suggests  that  we  can  obtain  a  filtered  version  of  the  WVD  by 
Fourier  transforming  the  Wxx(t;  s,  £). 
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IV.         FREQUENCY  HOPPED  SIGNALS  AND 
THEIR  CORRELATION  FUNCTIONS 

Communication  systems  can  utilize  a  large  number  of  digital  modulation  tech- 
niques. Among  those,  spread  spectrum  modulation  is  widely  used.  Spread  spectrum 
refers  to  any  modulation  scheme  that  produces  a  transmitted  bandwidth  much  larger 
than  the  information  bandwidth  [Ref.  50].  We  will  briefly  address  the  different  digital 
modulation  schemes  and  focus  on  frequency  hopping. 

A.      DIGITAL  MODULATION  SCHEMES 

For  comparison,  digital  communication  schemes  are  briefly  presented.  Digital 
modulation  techniques  use  the  binary  and  the  M-ary  schemes.  Binary  schemes  consist 
of  the  on-off  keying  (OOK),  also  called  ASK,  binary  phase  shift  keying  (BPSK),  and 
frequency  shift  keying  (FSK).  M-ary  schemes  are  generalizations  of  the  binary  schemes 
for  transmitting  M  symbols,  e.g.,  M-ary  PSK  or  M-ary  FSK. 

1.       Binary  schemes 

The  major  reference  of  this  topic  is  [Ref.  51]. 

•  OOK:  The  OOK  signal  is  represented  by: 

s(t)  =  Ac  m(t)  cos  wct,  (rV-1) 

where  Ac  is  the  carrier  signal  amplitude,  m(t)  is  the  transmitted  binary  data 
over  the  bit  duration  T&.  For  unipolar  representation  it  has  either  "1"  or  "0". 
Therefore,  over  an  observation  time  2],  the  signal  is  zero  (i.e.,  m(t)  =  0)  or  it 
is  a  single  sinusoid  at  frequency  uc  (i.e.,  m{i)  =  1). 

•  BPSK:  The  BPSK  signal  is  represented  by: 

s{t)  =  Ac  cos  [wct  +  A0m(t)} ,  (IV.2) 

where  m(t)  is  the  bi-polar  message  signal.  Over  the  bit  interval  Tb  the  signal 
has  the  value  ±1  depending  on  whether  the  bit  information  is  1  or  0.   Ad  is 
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the  corresponding  phase  shift.  Usually,  AO  is  chosen  as  tt/2.  Therefore,  the 
BPSK  signal  has  either  one  of  the  two  forms: 

s(t)  =  ±Acsinwct  ;      for    m(t)  =  ±1 .  (IV.3) 

Therefore,  under  any  value  of  m(i)  the  BPSK  signal  will  be  a  sinusoid  with  a 
fixed  frequency  component  but  with  a  phase  of  zero  or  7r. 

•  FSK:  The  continuous  phase  FSK  signal  is  given  by: 

s(t)  =  Ac  cos  [wct  +  6(t)] ,  (IY.4) 

and 

0(t)  =  Df[     m(£)di,  (TV. 5) 

where  Df  is  a  modulation  index.  For  binary  messages  m(t)  the  resultant 
FSK  signal  is  called  the  binary  FSK  scheme.  The  FSK  signal  has  one  of  two 
frequencies  without  phase  discontinuities  at  the  transition  points. 

2.       M-ary  Schemes 

For  an  M-ary  schemes  a  message  m(t)  has  M  symbols.  Consequently,  the 
transmitted  signal  will  have  M  different  states. 

•  M-ary  ASK:  The  M-ary  scheme  of  the  OOK  (or  ASK)  may  be  implemented 
for  different  values  of  M.  An  example  is  QAM  (i.e.,  M  =  4).  All  of  the  QAM 
states  have  the  same  single  frequency  components  but  differ  in  amplitude. 

•  M-ary  PSK  (MPSK):  The  M-ary  PSK  is  generated  similar  to  BPSK  with  the 
exception  that  the  value  of  AO  is  chosen  according  to  the  number  M.  Thus, 
all  the  MPSK  states  have  the  same  frequency  component  but  differ  in  phase. 

•  M-ary  FSK  (MFSK):  The  M-ary  scheme  is  similar  to  BFSK.  But,  instead  of 
two  symbols  (states),  it  has  a  set  of  M  symbols.  Consequently,  M  different 
frequencies  are  transmitted.  Thus,  the  MFSK  signal  can  be  expressed  as  s(t)  = 
Accos(wit  +  $i)  where  u>i  and  Oj  are  the  ith  frequency  and  phase  corresponding 
to  the  ith  symbol  of  the  message  m(t). 

In  summary,  digital  modulation  techniques  are  characterized  by  their  signaling 
mode  which  consists  of  a  single  fixed  frequency  component  over  the  duration  of  a  given 
bit. 
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B.      SPREAD  SPECTRUM  COMMUNICATION  SIGNALS 

The  major  reference  for  this  section  is  [Ref.  19,  50].  Spread  spectrum  (SS) 
communication  signals  are  characterized  by  a  wide  transmission  bandwidth,  and  a 
low  power  spectral  density.  SS  signals  have  two  main  advantages: 

1)  The  message  has  a  low  probability  of  being  intercepted  (LPI)  as  a  result  of 
the  wide  frequency  band,  and  the  low  power  spectral  density  of  the  signals. 

2)  SS  systems  can  reject  jamming  signals  and  allow  many  users  to  share  the  same 
frequency  band  as  a  result  of  the  spreading  gain  [Ref.  19,  50]. 

Among  the  different  possible  SS  modulation  formats,  the  following  are  preva- 
lent: 

1)  Frequency  Hopping  (FH):  The  complex  baseband  signal,  c(t),  with  basic  pulse 
shape  p(£),is  given  by 

c(t)  =  J2 exp  {j{2irfnt  +  <f>n] p[t  -  nTh],  (IV.6) 

where  Th  is  the  pulse  duration,  better  known  as  the  hop  interval.  The  pseudo- 
randomly  generated  sequence  of  frequencies  {fn}  will  drive  the  modulator  to 
generate  a  modulated  version  of  p(t).  {(f>n}  is  an  associated  phase  shift  at  each 
carrier  frequency  fn. 

2)  Direct  Sequence  (DS)  Modulation:  The  complex  baseband  signal,  c(t),  of  DS 
is  given  by 

c{t)  =  Y,CnP{t-nTc),  (IV.7) 

n 

where  {cn}  is  a  pseudorandom  sequence  which  modulates  a  sequence  of  pulses 
over  a  duration  Tc,  known  as  the  chip  interval. 

3)  Time  Hopping  (TH):  The  pulse  waveform  is  given  a  fraction  of  duration  T5, 
i.e.,  Ts/M.  A  typical  time  hopping  waveform  might  be 

eM«X>(«-(»  +  £)r.),  (IV.8) 

which  means  the  pulse  will  be  controlled  by  the  pseudo-random  number  (a„) 
to  appear  at  any  of  the  M  time  segments  within  the  duration  Ts. 

4)  Hybrid  Modulations:  A  blend  of  the  above  techniques  may  satisfy  better  per- 
formance according  to  some  requirements. 
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Figure  5  shows  a  typical  SS  system.  The  SS  spectrum  uses  two  types  of  mod- 
ulation. The  first  type,  the  baseband  data  modulation,  is  also  called  the  primary 
modulation.  The  second  type,  the  actual  SS  modulation,  is  also  called  the  secondary 
modulation.  For  simple  SS  system  realizations,  certain  combinations  of  primary  and 
secondary  modulations  are  usually  employed.  For  DS,  the  binary  phase  shift  keying 
(BPSK)  is  used  as  the  primary  data  modulation  scheme.  For  FH,  the  M-ary  frequency 
shift  keying  (MFSK)  is  used.  The  FH  system  which  implements  MFSK  data  modula- 
tion is  known  as  the  pure  FH  scheme,  and  is  the  most  popular  and  widely  applied  FH 
scheme.  Another  advantage  for  the  pure  FH  scheme  is  its  resistance  against  Repeat- 
back  jammers  (RBJ)  which  estimate  the  FH  signal  frequency  and  consequently  send 
a  proper  jamming  waveform.  Thus,  pure  FH  schemes  minimize  the  hostile  activity  of 
RBJ  as  each  hop  frequency  carries  one  symbol  of  the  transmitted  message  only. 

C.      THE  FREQUENCY  HOPPED  SIGNAL 

The  FH  of  the  BFSK  signal  is  given  by 

x(t)  =  V2Psin  [(uo  +  un  +  dnAcu)t] ,  (IV.9) 

where  P  is  the  signal  power.  The  frequencies  u>o  and  Au;  are  constants,  and  un  is  the 
frequency  for  the  nth  symbol  dn  whose  value  is  ±1. 

During  any  hop,  only  one  frequency  component  will  be  at  the  transmitter 
output  when  the  hopping  interval  equals  the  symbol  interval.  The  range  of  variation 
of  ljc  =  (cj«  +  u)q  ±  Aw)  is  known  as  the  hopping  bandwidth.  Two  types  of  FH 
systems  exist  depending  on  the  hop  interval  Th  and  the  symbol  interval  Ts:  fast- 
frequency  hopping  (FFH)  and  slow  frequency  hopping  (SFH)  which  differ  in  their 
hopping  speed.  The  number  of  hops  within  one  symbol  duration  is  the  measure  of 
the  speed  of  the  hopping  rate.  For  pure  FH  signals  with  BFSK,  the  SFH  requires 
Th  >  Ts  while  FFH  requires  Th<Ts.  The  SFH  contains  one  or  more  symbols  during 
the  hop  interval  while  the  FFH  has  one  or  more  hops  over  one  symbol  interval.  Figure 
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6  shows  a  typical  pure  FH  system  using  the  BFSK  as  a  primary  modulation.  Figure 

7  illustrates  a  typical  time-frequency  behavior  of  the  BFSK  pure  FH  scheme. 

D.      THE  INSTANTANEOUS  CORRELATION  FUNCTION 
OF  FREQUENCY  HOPPED  SIGNALS 

Spread  spectrum  studies  have  usually  considered  the  FH  signal  as  a  stationary 
process  [Ref.  19,  50]  even  though  the  spectrum  of  the  FH  signal  varies  with  each 
hop  interval.  Therefore,  the  correlation  representation,  using  time  averaging,  is  not 
suitable  for  this  nonstationary  process. 

One  way  to  identify  the  FH  signal  is  to  monitor  the  time  evolution  of  the  signal. 
Hence,  we  need  to  keep  the  time  dependency  in  the  correlation  representation.  In 
this  work  we  resort  to  the  time-varying  correlation  definition  for  characterizing  the 
FH  signals.  Therefore,  we  select  the  instantaneous  correlation  function  (ICF)  as  the 
candidate  correlation  representation.  In  this  section,  we  will  address  the  structure  of 
the  ICF  of  the  pure  FH  signals. 

The  pure  FH  signal  may  be  represented  as  successive  intervals  (i.e.,  hops) 
with  single  complex  exponentials.  The  frequency  within  each  interval  (i.e.,  the  hop 
frequency)  is  controlled  by  a  randomly  generated  sequence  of  numbers,  therefore,  we 
can  assume  the  following  about  the  used  frequencies: 

(1)  The  number  of  the  different  hopping  frequencies  is  much  larger  than  the  num- 
ber of  observed  hops. 

(2)  The  selection  of  discrete  hopping  frequencies  (hopping  pattern)  is  determined 
by  a  pseudo-random  number  generator  (PNG),  i.e.,  all  hopping  frequencies  are 
equi-probable. 

(3)  The  observation  interval  is  much  smaller  than  the  period  of  PNG.  The  period 
of  PNG  is  the  number  of  generated  hops  times  the  hop  interval. 

As  a  result,  we  can  deduce  that  any  two  successive  hops  will  have  two  dif- 
ferent frequencies.    The  difference  in  the  frequency  of  adjacent  hops  will  generate 
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the  structure  of  the  instantaneous  correlation  functions.  Recall  that  we  define  the 
instantaneous  correlation  function  as: 

where  t  is  time  and  r  is  time  lag.  Over  the  two-dimensional  plane  of  time  and  time 
lag  (i.e.,  t  —  t  plane)  we  compute  the  ICF  for  the  time  lag  \r\  <  7^,  which  will  allow 
correlating  adjacent  hops  only.  If  the  values  of  (t  +  |)  and  (t  —  |)  are  both  confined 
within  the  Lth  hop  then  the  ICF  will  be  given  by: 

R{t,r)    =    eJ"L(t+T/2)  ^ut-r/vy 

=    ejuJLT,  (IV.  11) 

which  is  a  function  of  r  and  ljl  only.  Note  that  the  values  of  (t  +  |)  and  (t  —  | )  are 
both  confined  within  the  same  Lth  hop  if  they  satisfy, 

(L  -  \)Th  <{t  +  \)  and  (t  - 1)  <  LTh.  (IV.12) 

This  inequality  forms  the  boundaries  of  the  diamond  cellular  shape  for  different 
values  of  L.  This  cellular  structure  is  shown  in  Figure  8.  Inside  each  diamond  DL, 
the  ICF  results  from  correlating  signals  belonging  to  the  same  hop,  while  outside  the 
diamond  the  ICF  results  from  correlating  signals  belonging  to  two  consecutive  hops. 

A  detailed  derivation  of  the  ICF  is  given  in  Appendix  A,  which  leads  to  the 
doubly  indexed  correlation  function  Rm,n(t,T)  given  by: 

-Rm,n(*,  r)  =  exp  j  I  (um  -  un)t  +  (um  +  un)-  J  .  (IV.13) 

where  m  and  n  are  the  indices  of  the  two  adjacent  hops.  Note  that: 

1)  Within  the  main  diamond  of  the  nth  hop,  i.e.,  m  =  n,  the  ICF  is  given  by 

Rn,n(tiT)  =  exp{i(o;nr)}.  (IV.14) 
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2)  Outside  the  main  diamond  in  the  upper  triangle  between  hop  numbers  (n)  and 
(n  +  1),  the  ICF  is  given  by 

R„+i,„(t,  t)  =  exp  |j(u;n+1  -  un)t  +  {un+i  +  u)n) -  J  .  (IV.15) 

In  conclusion,  the  instantaneous  correlation  function  of  a  pure  FH  signal  has 
a  cellular  or  diamond  structure.  Inside  the  Lth  diamond  it  has  a  single  complex  expo- 
nential component  along  the  delay  axis  representing  the  Lth  hop  frequency.  Outside 
the  diamond,  R(t,r)  is  a  product  of  two  terms,  e^Un~Um>it  and  eJ(Wn+Wm)T/2,  where  ujn 
and  um  are  two  consecutive  hop  frequencies. 
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Figure  7.  Time-Frequency  behavior  of  BFSK-FH  signal. 
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Figure  9.  Areas  and  indices  of  the  doubly  indexed  function  £,„,„(*,  r). 
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V.         ANALYSIS  USING  MORLET  WAVELETS 

We  studied  the  structure  of  the  ICF  surface  obtained  for  complex  FH  signal 
in  chapter  IV,  and  showed  that  it  consists  of  complex  sinusoidal  components.  In  this 
chapter,  we  will  analyze  the  Morlet  wavelet  transform  of  the  ICF.  We  will  use  the 
Morlet  basis  function  for  two  reasons.  First,  the  Morlet  wavelet  is  a  complex  sinusoid, 
modulated  by  a  Gaussian  window,  and  is  inherently  best  suited  for  filtering  sinusoidal 
signals.  Second,  the  Morlet  wavelet  has  a  mathematical  formulation  which  makes  the 
analysis  and  derivations  tractable. 

Recall  that  the  mother  Morlet  wavelet  g(t)  is  given  by: 

g(t)  =  ejkte-t2/2a\     for  -  oo  <  t  <  +00 ,  (V.l) 

where  A:  is  a  constant  that  represents  the  modulation  frequency,  and  a2  is  inversely 
proportional  to  the  roll-off  factor  of  the  Gaussian  window.  The  Fourier  transform  is 
given  by: 

G{u)  =  e~^-k)2 .  (V.2) 

The  Morlet  wavelet  does  not  satisfy  the  admissibility  nor  the  orthogonality 
condition  as  it  exists  over  an  infinite  time  interval,  but  for  practical  applications  it  is 
truncated  when  it  is  sufficiently  attenuated.  If  the  wavelet  has  finite  support  (non-zero 
only  over  finite  interval),  then  the  scaled  and  shifted  wavelet,  gs,i{t),  will  be  non-zero 
over  the  interval  [£  —  nas,£  +  nas],  where  s  is  the  scale  and  I  is  the  shift.  Here,  n 
is  a  preselected  number  which  ensures  sufficient  decay.  The  functional  dependence  of 
gs,e(t)  is  as  described  in  Equation  11.41.  Usually,  n  is  chosen  >  4  [Ref.  46]. 

A.      TRANSFORM  OF  A  TIME-LIMITED  COMPLEX 
EXPONENTIAL 

We  showed  in  Chapter  IV  that  the  ICF  of  the  FH  signal  consists  of  complex 

exponential  components  with  frequencies  that  differ  from  region  to  region.  Therefore, 
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we  need  to  address  the  response  of  the  wavelet  transform  to  complex  exponential 
signals.  Let  the  signal  f(t)  be  a  complex  exponential  Ae^Lt  defined  over  the  interval 
t  £  [LT,  (L  +  1)T],  where  T  is  the  hop  interval,  and  L  stands  for  the  Lth  hop.  When 
<7  =  1,  the  Morlet  wavelet  transform  (MWT)  is  given  by: 

1        ri+ns  o 

Wis,  £)  =  a4=  e3«t.te-JW-me-W-tf/*)dt  t  (v.3) 

V  s  Je-ns 

where  £  is  confined  to  the  interval  [LT  +  ns,  (L  +  V)T  —  ns\.  This  interval  guarantees 
that  the  wavelet  will  be  confined  over  the  same  hop  without  crossing  the  border  to 
another  hop.  Later  on  we  will  compute  the  wavelet  transform  when  the  signal  runs 
over  the  border  between  two  adjacent  hops.  The  exponent  in  Equation  V.3  can  be 
written  as: 


It 


£ 


k£    \e 


-«T^+  [3VL-3-  +  -Z  )t  +  j 


2  5 

Substituting  coL  —  k/s  =  u  we  obtain 

W(sJ)  =  ,4-I=e-('2/2s2)+jWs)  r" 

y/S  Je-ns 

This  equation  can  be  reduced  to 


exp 


2s1 


2s< 


t2+[ju  +  -)t 


dt 


(V.4) 


(V.5) 


W(s,£)  =  A]f^exp 


52w£      k 


~~2+  Jw^  +  iw^ 


(EF), 


(V.6) 


where  EF  is  given  by: 


EF  =  erf 


-^={n  +  j{uLs-k)) 


+  erf 


V2 


(n-j(ujLs-k)) 


(V.7) 


for  £e[LT-r  sn,  (L  +  1)T  -  sn]. 

We  note  that  the  magnitude  of  W(s,  £)  is  independent  of  the  wavelet  shift 
variable  £. 

The  term  EF  defined  above  can  be  approximated  using  the  following  formula  for  the 
complex  error  function  [Ref.  54]  with 
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-X2 

e  x 


erf(x  +  jy)    =    erf(x)  +  - — [1  —  cos  2xy  +  j  sin  2xy] 

Z7TX 

2  _  2   °°    e~(1/4)m2 
+"e  X   £  m2  +  4:r2  [/m(g»  ri  +  J9m{x,  y)] 

+e(x,y),  (V.8) 

where 

/m(#j  y)    =    2a:  —  2a:  cosh  (my)  cos(2a:y)  +  m  sinh  (my)  sin(2a:y) 
9m{x,  y)    —    2a:  cosh(ray)  sin(2a:y)  +  m  sinh(my)  cos(2a:y)  , 

and  \e(x,y)\  «  10"16|er/(a:,  y)|. 
Since  .EF  is  of  the  form 

£F(a:,  y)  =  er/(a:  +  jy)  +  er/(x  -  jy), 

we  need  to  determine  erf(x  —  jy).  Recall  that  sin(x)  and  sinh(a:)  are  odd  functions 
while  cos  (a:)  and  cosh  (a:)  are  even  functions.  Thus  we  can  express 


-X2 

e  x 


erf(x  -  jy)    =    erf(x)  +  - —  [(1  -  cos  2xy)  -  j  sin  2a:y] 


2  _  2  °°    e~(1/4)m* 

+~e  ■  £  m2  +  4;r2  (/»(«» -y) + ifl»(*»  -y)) » 

where 

fm{x,—y)    =    2x  —  2x  cosh(my)  cos(2a:y)  +  msinh(my)sin(2a:y), 
<?m(a:,  — y)    =    —  2a:  cosh  (my)  sin(2a:y)  —  msinh(my)  cos(2a:y) . 

Or,  equivalently, 

fm(x,-y)   =   fm(x,y), 

9m{x,-y)      =      -0m(*,2/). 
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(V.9) 


Using 


and 


erf(x  -  jy)  =  erf*{x  +  jy), 


erf(x  +  jy)  +  erf(x  -  jy)  =  2Re  {erf(x  +  jy)} , 


(V.10) 


(V.H) 


where  Re{-}  denotes  the  real  part,  we  can  express  W(s,£)  as: 


W(sJ)    =    AV2ns  exp 


2,  ,2 


S  CO 


-  —  +  su)Lk  +  juL£ 


1 


Re  <  erf  (  -j=[n  +  j(o;Ls  -  A)]  J  \  , 


(V.12) 


for  /  €  [LT  +  ns,  (L  +  1)T  -  ns\. 

Equation  V.12  shows  that  the  transform  has  a  linear  phase  response,  $(a;), 
independent  of  the  scale  s,  which  is  given  by: 


$(w)  =  ujl£. 


(V.13) 


Note  that  if  the  complex  exponential  signal  has  an  initial  phase  shift,  $l,  the  phase 
of  the  output  of  the  transform  will  change  to 


$(u)  =  {uL£  +  0L). 


(V-14) 


The  amplitude  may  be  expressed  as  AB(s,uj),  where  B(s,u))  is  the  wavelet 
gain  given  by: 


i?(s,a;)  =  v2tts  exp 


s2u2L      k2 

—  -1  +  suLk 


J^lerf{-^=[n  +  j{coLs-k)]j\. 


(V.15) 

The  gain  of  the  Morlet  wavelet  transform  B(s,co)  is  plotted  in  Figure  10  as  a  function 
of  scale  and  frequency. 

The  spectral  response  of  the  wavelet  is  plotted  for  values  of  s  =  0.5,1,  and 
2.  Over  this  dyadic  grid  the  spectral  response  contains  some  regions  of  low  response 
due  to  the  narrow  bandwidth  of  the  Morlet  wavelet  filters.  In  many  applications, 
discrete-scale  wavelets  are  preferred.  Note  that,  the  scale  grid  needs  to  be  sampled 
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linearly  (i.e.,  s  =  1, 2,3,4,  ... )  or  the  roll-off  factor  (cr2)  of  Equation  V.l  has  to  be 
decreased  when  applying  Morlet  wavelets  for  full  spectral  coverage. 
Complex  error  function.     The  error  function  used  in  W(s,£)  has  a  complex  ar- 
gument and  we  use  a  computational  method  given  in  [Ref.  54,  56]  to  compute  it. 
The  error  function  is  defined  as 

erf(z)  =-=  [*  e~u2du,  (V.16) 

y/TT  JO 

while  the  complementary  error  function  is  defined  as: 

erfc(z)  =  1  -  erf(z).  (V.17) 

Let  w(z)  be  defined  as: 

w(z)  =  e~z2erfc(-jz) .  (V.18) 

Then,  for  a  complex  argument  z,  the  error  function  is  expressed  as: 

erf(z)  =  1  -  e~z2w(jz).  (V.19) 

The  function  w(z)  is  tabulated  for  some  values  of  z  in  [Ref.  54]  and  may  be 
computed  using  the  algorithm  described  in  [Ref.  56]. 

B.      ANALYSIS  OF  THE  TRANSITION  REGION 

The  transition  region  is  defined  as  that  where  the  FH  signal  f(t)  changes  from 
one  to  another  hop  at  time  tr.  Thus,  the  wavelet  transform  involves  two  different 
signal  frequency  components  in  the  transition  region.  The  transition  region  between 
interval  L  and  interval  L  +  1  will  be  defined  as  £  €  \{L  +  1)T  —  ns,  (L  +  1)T  +  ns], 
where  the  Morlet  wavelet  support  interval  is  2ns. 

Hence,  the  Morlet  wavelet  transform  is  given  by: 

where  tr  —  ns  <  £  <  tT  +  ns. 
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The  two  parts  of  Equation  V.20  can  be  solved  similarly  to  Equation  V.3,  which 


leads  to 


W(s,£)    =    y^exp 
{erf 


S2"l        ^  ,         •       n 

— 7T-  -  -TT  +  SULK  +  JU)L£ 


V2 


-KS 

+VTexp 


J{uLs-k) 


+  erf  (  —j=[n  +  j(uLs  -  k)] 


s2ul+1      k2  . 

— Y~  -  y  +  suL+lk  +  jWL+l« 


er/(-^[n-j(o;L+i5-fe)] 


erf  17! 


6f  C  . 


-j{cuL+iS-k) 


(V.21) 


Over  the  length  of  the  transition  region  the  magnitude  of  MWT  depends  on 
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Figures  11  and  12  show  an  example  of  the  magnitude  and  the  phase  of  the  wavelet 
transform  over  a  transition  region.  Two  frequencies  are  used  (4  and  7  rad/sec)  and 
the  magnitude  and  the  phase  of  the  wavelet  transform  are  plotted  for  scale  8  =  1. 
The  magnitude  of  the  wavelet  coefficients,  for  the  signal  with  the  first  frequency  (i.e., 
f=4) ,  is  larger  than  the  magnitude  of  the  wavelet  coefficients  of  the  signal  with  the 
second  frequency  (i.e.,  f=7).  The  magnitude  and  the  phase  of  the  wavelet  transform 
for  scale  s  =  0.5  is  plotted  in  Figure  12.  The  magnitude  of  the  wavelet  coefficients, 
for  the  signal  with  the  first  frequency  (i.e.,  f=4),  is  smaller  than  the  magnitude  of 
the  wavelet  coefficients  of  the  signal  with  the  second  frequency  (i.e.,  f=7).  Note  also 
that  there  is  a  gradual  transition  between  the  two  magnitudes  according  to  Equation 
V.21.  This  indicates  that  the  Morlet  wavelet  transform  does  not  produce  spikes  at 
the  transition  points  between  frequency  hops. 
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C.      ANALYSIS  OF  IDEAL  WHITE  GAUSSIAN  NOISE 

We  investigate  the  response  of  the  Morlet  wavelet  to  ideal  white  Gaussian 
noise.  The  WT  is  given  by: 

Wn(s, I)  =  r  l-n(t)g*(*—*)dt.  (V.22) 

J  —  oo  y/S  S 

For  white  Gaussian  noise  n(t)  and  due  to  the  linearity  of  the  wavelet  transform,  the 
output  Wn(s,  £)  has  a  Gaussian  distribution.  Since  n{i)  is  assumed  having  zero  mean 
and  variance  o\,  then  the  mean  value  of  Wn(s,£)  is 

E{Wn(sJ)}    =    sf^intfjc-^^^e-W^^j 

=    0.  (V.23) 

Since  Wn(s,£)  is  zero  mean,  its  variance  is  given  by: 
°l    =    E{\W(sJ)\2}  =  E{W(S,£)W*(SJ)} 

1     f°°      f°°     ,       ,      m  fr     /      x     */      m      ftl-t\      *{h-2\] 


S  J—oo  J—oo  \        S        I 


(V.24) 

When  two  samples  n{t\)  and  nfa)  are  assumed  independent  identically  distributed 
Gaussian  (i.i.d.)  random  variables,  the  term  E{[n(ti)n*(t2)}  is  equal  to  zero  except 
when  ti  =  t2.  Hence  we  have 

r2  ^00 


2       ^  r 

b    J  —c 


■t-t    2 


dt 


=    a2nV^.  (V.25) 

Consequently,  the  Morlet  wavelet  transform  of  white  Gaussian  noise,  with  zero 
mean  and  a\  variance,  is  Gaussian  with  zero  mean  and  a  variance  of  (J^V^  at  any 
wavelet  scale  s.  Note  that  we  assumed  an  infinite  support  for  the  Morlet  wavelet  in  the 
above  derivation.  If  we  carry  out  the  same  procedure  for  the  finite  support  wavelet, 
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we  will  have  a  different  expression  for  the  variance  of  Wn(s,£)  due  to  integrating  the 

ry 

term  \g  (— )     over  a  finite  interval.  In  such  a  case,  the  variance  becomes 


2         °n    I* 

S     J  —7 


't-t    2 


dt,  (V.26) 


which  results  in 

a2w(sj)=a2n^erf(ns).  (V.27) 

We  note  that  for  ns  »  1,  erf(ns)  is  well  approximated  by  1.  Thus,  we  conclude  that 
the  variance  of  the  wavelet  transform  will  be  essentially  independent  of  the  wavelet 
scale  and  can  be  approximated  by  o^y/ft. 

Actual  Noise  of  the  Wavelet  Surfaces.  Application  of  the  wavelet  transform 
to  the  ICF  surface  will  result  in  a  number  of  wavelet  surfaces  corresponding  to  the 
number  of  wavelet  scales  used.  Since,  we  exploit  the  wavelet  surfaces  for  identify- 
ing the  FH  signal,  we  need  to  investigate  the  noise  distribution  over  these  surfaces. 
The  noise  of  the  wavelet  surfaces  is  approximated  as  Gaussian  noise.  The  following 
considerations  assist  in  making  this  decision: 

1)  The  noise  background  is  assumed  to  be  additive  white  Gaussian  noise.  The  FH 
sinusoidal  components  and  the  noise  are  assumed  to  be  independent.  Noise 
realizations  at  two  different  time  instants  are  assumed  to  be  uncorrelated  and 
independent. 

2)  The  ICF  surface  consists  of  noise  components  which  result  from  the  product  of 
two  independent  Gaussian  random  variables.  This  resulting  noise  surface  has  a 
A'-distribution  shape  shown  in  Figure  13,  where  r  is  the  correlation  coefficient 
of  the  two  Gaussian  random  variables  [Ref.  57]. 

3)  The  wavelet  transform  of  the  ICF  surface  is  a  weighted  sum  of  the  ICF  samples. 

According  to  the  central  limit  theory,  usually,  the  sum  of  number  of  identically 
distributed  random  variables  tends  to  a  Normal  (Gaussian)  distribution.  Since  we 
take  the  wavelet  transform  of  the  ICF  surface  in  the  time  delay  direction,  the  samples 
of  the  wavelet  surfaces,  in  the  time  delay  direction,  may  be  approximated  as  having 
a  Gaussian  distribution.    Simulation  results  show  that  the  noise  distribution  over 
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the  wavelet  surfaces,  taken  in  direction  of  time  delay,  has  approximately  a  Gaussian 
distribution.  The  Gaussian  assumption  was  tested  using  a  chi-square  test.  The  main 
central  part  of  the  noise  distribution  (excluding  the  tails)  was  found  to  be  Gaussian 
with  confidence  level  >  90%. 
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Figure  10.  Morlet  wavelet  gain. 
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Figure  11.  Morlet  wavelet  gain  and  phase  over  the  transition  region  (s=l). 
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Figure  12.  Morlet  wavelet  gain  and  phase  over  the  transition  region  (s=0.5). 
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Figure  13.  The  K-distribution. 
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VI.         ANALYSIS  OF  THE  PROCESSING 

SCHEME 

A.      INTRODUCTION 

Our  goal  is  to  exploit  the  wavelet  transform  of  the  ICF  surface  in  order  to 
identify  FH  modulation  schemes.  The  wavelet  transform  results  in  a  set  of  wavelet 
surfaces,  one  for  each  scale.  We  address  the  interception  problem  using  two  ap- 
proaches, either  of  them  can  be  used  independently,  but  they  differ  in  performance  as 
will  be  shown  later.  In  the  first  approach,  we  visually  inspect  the  2-D  wavelet  surfaces 
to  identify  and  classify  the  structure  of  the  FH  signal  and  obtain  an  estimate  for  the 
hop  time  interval.  In  the  second  approach,  we  apply  a  proposed  processing  scheme  to 
estimate  the  hop  start/stop  times,  the  hop-scale  pattern,  and  the  hop  frequency.  The 
extraction  of  the  hop  start/stop  times  is  addressed  using  an  edge  detection  approach 
by  applying  a  compass  operator  which  is  well  known  in  the  image  processing  area. 
The  hop-scale  pattern  is  obtained  by  applying  an  energy  analysis.  The  energy  anal- 
ysis assigns  a  scale  index  (called  the  proper  scale)  to  each  hop.  The  proper  scale,  for 
each  hop,  is  that  scale  which  has  the  greatest  energy  content.  The  sequence  of  proper 
scales,  representing  the  hop  sequence,  is  called  the  hop-scale  pattern.  The  frequency 
of  each  hop  can  be  extracted  from  the  wavelet  surface  at  the  proper  scale. 

The  identification  of  the  FH  signal  may  be  accomplished  based  on: 

(1)  The  hop-scale  pattern:  If  two  or  more  wavelet  scales  are  applied,  the  hop-scale 
pattern  of  the  FH  signal  will  be  different  from  the  hop-scale  patterns  of  other 
digital  modulation  signals. 

(2)  The  frequency  diversity:  If  all  frequencies  reside  in  one  scale,  then,  an  FH  sig- 
nal will  have  different  frequency  components  as  a  function  of  the  hop  intervals. 

Recall  that,  the  ICF  of  the  FH  signal  has  a  cellular  structure  which  comprises 
of  a  tiling  of  diamonds.  Each  hop  results  in  a  diamond  with  a  width  equal  to  the  hop 
interval.  Thus,  the  diamond  boundaries  point  to  the  hop  start/stop  times.  All  wavelet 
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surfaces  at  all  scales  emphasize  the  signals  and  hence  the  diamonds  which  belong  to 
that  scale.  Other  signals  are  attenuated.  In  this  chapter,  we  address  aspects  of  the 
discrete  implementation,  the  proposed  processing  scheme,  the  measured  parameters, 
and  performance  measures.  In  section  VLB,  we  address  the  discrete  form  of  the 
ICF  and  its  output  expression  to  a  white  noise  input.  In  section  VI. C,  we  address 
the  visual  inspection  and  FH  identification  from  the  wavelet  surfaces  and  obtain  an 
estimate  for  hop  time  intervals.  Then  we  compare  the  FH  wavelet  surfaces  with 
wavelet  surfaces  from  other  digital  modulation  schemes.  In  section  VI.D,  we  consider 
an  energy  analysis  for  identifying  the  scale  of  each  frequency  hop.  We  also  investigate 
the  performance  of  scale  identification  and  evaluation  measures.  In  section  VI.E, 
we  address  the  equalization  of  the  spectral  shape  of  the  ICF  and  its  impacts  on  the 
performance  of  scale  identification.  In  section  VI.F,  we  address  the  hop  frequency 
extraction  from  the  wavelet  surfaces  at  the  proper  scale  and  from  the  original  time 
signal.  We  also  investigate  the  performance  of  frequency  extraction  and  evaluation 
measures.  Finally,  we  address  the  extraction  of  the  hop  start/stop  times  in  section 
VI.G. 

B.      PROCESSING  SCHEME 

The  interception  problem  usually  assumes  some  prior  knowledge  about  the  sig- 
nal of  interest.  In  our  case,  we  assume  the  signal  hopping  bandwidth  is  approximately 
known  and  the  data  is  properly  heterodyned  and  sampled. 

1.       Processing  Sequence 

The  input  data  parameters  to  our  processing  scheme  are  the  parameters  of  the 
FH  signal  (i.e.,  sampling  frequency  and  the  range  of  possible  frequencies)  while  the 
outputs  are: 

(a)  The  wavelet  transform  of  the  ICF  surface,  i.e.,  the  wavelet  surfaces. 

(b)  The  scale  of  each  hop. 

(c)  The  frequency  of  each  hop. 
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(d)  The  hop  interval  start/stop  times. 

The  data  (input)  process  generates  FH  signals  according  to  some  predeter- 
mined input  parameters  (i.e.,  number  of  hops,  hop  interval,  hop  frequencies,  sampling 
frequency  and  SNR).  The  analysis  process  generates  the  ICF  surface  and  computes 
its  wavelet  transform  at  the  predetermined  wavelet  scales.  The  measurement  process 
extracts  the  outputs  (the  scale  of  each  hop,  the  frequency  of  each  hop,  and  the  hop 
interval  start/stop  times). 

Figure  14  shows  the  functional  description  of  the  processing  scheme  where  the 
Input  stands  for  the  input  process.  The  Hilbert  Transform  generates  the  analytic  form 
of  the  signal.  The  Instantaneous  Correlation  Function  generates  the  (discrete  time) 
ICF  surface.  The  Wavelet  Transform  computes  the  (discrete  time)  wavelet  transform 
of  the  ICF  surface  (i.e.,  wavelet  surfaces).  The  Scale  Identifier  identifies  the  hop  scale. 
The  Hop  Frequency  extracts  the  hop  frequency  from  the  wavelet  surfaces.  The  Hop 
Timing  extracts  the  hop  start/stop  times  and  the  hop  interval. 

2.       Discrete-Time  Implementation  of  the  Instantaneous 
Correlation  Function 

For  discrete  implementations,  let  R(n,  u)  define  the  ICF  of  the  discrete-time 

signal  x(n), 

R{n,  u)  =  x  (n  + 1)  x*  (n  -  |)  ,  (VI.l) 

where  n  is  time  and  u  is  the  time  delay.  Note  that,  the  normalized  sampling  interval 
is  1  for  the  discrete  implementation.  In  addition,  the  combined  index  n  ±  |  should 
assume  only  integer  values.  Therefore,  there  will  be  two  choices: 

•  Zero-inserted  ICF:  If  we  let  u  assume  only  even  integer  values,  then  R(n,u), 
for  odd  u,  must  be  set  to  zero. 

•  Frequency-doubled  ICF:  If  we  let  u  —  2m  then  R(n,  m)  will  be  defined  as 

R(n,m)  =  x{n  +  m)x*(n  —  m) .  (VI. 2) 
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Zero-inserted  ICF.  The  effect  of  zero  insertion  to  the  ICF  will  result  in  doubling 
the  number  of  wavelet  coefficient  over  the  resultant  surfaces,  implying  more  computa- 
tion and  storage  load.  Each  hop  will  reside  within  its  scale  according  to  its  associated 
signal  frequency. 

Frequency- doubled  ICF.  The  frequency-doubled  discrete  ICF  has  the  following 
characteristics  in  contrast  to  the  zero-inserted  ICF: 

•  Due  to  doubling  the  frequency,  the  minimum  sampling  frequency  of  the  mon- 
itored signal  will  be  twice  the  Nyquist  rate. 

•  We  obtain  half  the  number  of  coefficients  at  the  ICF  surface  and  consequently 
at  the  wavelet  surfaces. 

•  Over  the  wavelet  surfaces,  each  hop  will  be  located  at  a  lower  scale  index  than 
its  associated  signal  frequency  (lower  scale  index  indicates  higher  frequency 
according  to  the  MATLAB  designation,  see  VI.A.3). 

Since  we  usually  pick  a  sampling  frequency  which  is  a  multiple  of  the  Nyquist 
rate,  the  choice  of  the  frequency  doubled  ICF  is  advantageous  because  of  a  fewer 
number  of  coefficients,  hence,  saving  in  computation  and  storage.  On  the  other  hand, 
the  loss  of  processing  gain,  due  to  pushing  each  hop  to  a  lower  scale  index,  can  be 
overcome  by  chosing  higher  sampling  frequency.  For  these  reasons  we  adopted  the 
definition  of  the  frequency-doubled  ICF  in  the  processing  scheme. 

3.       MATLAB  Discrete  Wavelet  Transform 

The  MATLAB  Wavelet  Toolbox  is  used  to  implement  the  discrete  wavelet 
analysis  of  the  ICF  surface.  The  wavelet  transformation  is  carried  out  using  the 
one-dimensional  wavelet  transform  in  the  direction  of  the  time  delay  u  for  each  time 
element  of  the  correlation  function  (as  previously  addressed  in  Chapter  III).  MATLAB 
uses  a  scale  index  designation  and  a  corresponding  multi-level  decomposition  tree 
as  shown  in  Figure  15.  Sj  denotes  the  output  of  the  wavelet  transform  at  the  ^"th 
wavelet  scale  (i.e.,  the  wavelet  coefficients  or  the  details  at  the  jth  level).  Aj  denotes 
the  approximation  at  the  jth  level.  H0  is  the  scaling  filter  (low  pass  filter)  and  Hi  is 
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the  wavelet  filter  (high  pass  filter)  [Ref.  59].  MATLAB  uses  the  Daub-N  designation 
where  N  indicates  the  wavelet  order  and  the  number  of  vanishing  moments,  while  the 
actual  filter  length  is  2N  [Ref.  59]. 

C.      VISUAL  IDENTIFICATION 

In  this  section,  we  address  the  visual  inspection  of  the  wavelet  surfaces  in 
order  to  identify  the  FH  modulation  scheme.  We  investigate  different  approaches  to 
represent  these  surfaces  and  considered  real  and  complex  valued  wavelet  functions. 

1.       Real  and  Complex- Valued  Wavelets 

According  to  Daubechies  [Ref.  42],  there  is  no  symmetric  or  anti-symmetric 
real  compactly-supported  orthonormal  wavelet.  Symmetry  property  can  be  achieved 
only  for  complex- valued  wavelets.  Symmetry  of  a  wavelet  implies  that  the  FIR  filter 
representation  has  a  linear  phase  response  [Ref.  36].  It  is  an  important  feature  in  some 
wavelet  applications,  such  as  numerical  resolution  of  partial  differential  equations 
with  boundary  conditions  [Ref.  60].  A  complex- valued  Daubechies  wavelet  of  order 
3  is  given  in  [Ref.  60]  and  has  been  used  for  comparison  with  the  MATLAB  real 
Daubechies  wavelet  of  the  same  order.  The  Daub-3  complex  wavelet  has  the  following 
scaling  coefficients: 


h0(l)  =  h0(6)    = 
ft0(2)  =  M5)    = 


3  +  M5 
64 

h-jy/TE 
64 


M3)  =  M4)    =     ljL±ifl-  (VI-3) 

2.       Surface  Representation 

Structure  of  the  FH  signal  can  be  identified  and  classified  by  visually  inspecting 
the  wavelet  surfaces.  We  can  also  obtain  an  estimate  for  the  hop  start /stop  times.  The 
next  automated  step  is  achieved  by  applying  a  proposed  processing  scheme  to  estimate 
the  hop  start/stop  times,  the  hop-scale  pattern,  and  the  hop  frequency.   Using  the 
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analytic  form  of  the  FH  signal,  the  ICF  surface  will  be  a  complex- valued  function. 
Therefore,  for  any  real  or  complex  wavelet  function  the  resultant  two-dimensional 
wavelet  surfaces  are  complex- valued.  Hence,  we  consider  the  real  part,  the  imaginary 
part,  the  magnitude,  or  the  phase  (angle)  of  the  complex- valued  surface. 

An  opinion  test  was  carried  out  to  assess  the  quality  of  identification  from  the 
wavelet  surfaces  using  different  representations.  The  Morlet  wavelet  and  Daubechies 
wavelet  of  order  3  were  both  used  in  their  real  and  complex  form.  Four  different 
representations  were  investigated:  real  part,  imaginary  part,  magnitude  (absolute 
value),  and  phase.  The  objective  of  the  opinion  test  is  to  identify  the  cellular  structure 
over  the  wavelet  surface  and  to  identify  the  diamond  boundaries  as  an  indication  of 
the  hop  start/stop  times. 

Figures  16  to  19  present  the  real  part  of  the  WT  obtained  with  the  real-valued 
Daubechies  wavelet  of  order  3,  at  an  input  SNR  values  of  10  and  3  dB.  In  Figures 
16  and  17  we  note  that  we  can  not  identify  the  FH  structure  over  the  surface  of  the 
ICF  denoted  by  "CF".  For  wavelet  surfaces,  labeled  "Sir",  "S2r",...,  "S5r",  we  can 
identify  diamond  patterns  at  the  hops  number  1,  2,...,  5  ,  respecively.  Same  findings 
can  be  observed  in  Figures  18  and  19  about  wavelet  surfaces.  These  results  show  that 
we  can  identify  the  FH  structure  from  the  wavelet  surfaces,  while  it  is  not  possible  to 
do  so  from  the  ICF  surface. 

D.      ENERGY  ANALYSIS  AND  SCALE 
IDENTIFICATION 

Energy  analysis  is  performed  over  the  wavelet  surfaces  at  all  scales  considered. 

For  the  energy  analysis  we  assume  correct  hop  timing.    The  hop  interval  and  hop 

start/stop  times  are  estimated  by  visual  inspection  (and  later  on  from  the  processing 

scheme).  Thus,  we  can  point  to  each  hop  over  the  wavelet  surface  and  compute  the 

energy  contained.    The  energy  contained  at  all  scales  are  compared  to  identify  the 

proper  scale  of  each  hop. 
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1.        Energy  and  Energy  per  Sample 

Parseval's  theorem  for  the  complete  orthogonal  filter  bank  (over  L  partitions) 
is  applied  to  the  discrete  time  wavelet  analysis  [Ref.  45].  It  is  given  by: 

IW»)II2=E  \\C(L,2k)\2  +  J:\d(j,2k  +  l)A,  (VIA) 

kez  \  j=l  J 

where,  in  the  sense  of  wavelet  analysis,  C(L,  2k)  are  the  scaling  coefficients  at  the 
scale  L,  d(j,  2k  +  1)  are  the  wavelet  coefficients  at  scale  j,  and  k  is  the  wavelet  shift 
variable. 

Therefore,  the  quantity 

oo 

m=  e  w.*)ia  (vl-5) 

k=— oo 

represents  the  portion  of  the  signal's  energy  over  the  jth  scale.  For  narrowband  signals 
which  reside  within  one  scale,  the  signal  energy  will  belong  to  that  scale.  However, 
some  portions  of  the  signal  energy  will  be  leaking  to  other  scales.  To  identify  the 
proper  scale  one  can  consider  the  maximum  value  due  to  the  following  quantities: 

1)  Wavelet  coefficient. 

2)  Total  energy. 

3)  Energy  per  sample. 


Energy  per  sample  is  defined  as: 


MS)  =  fg,  (VI'6> 


where  j  is  the  scale  index,  E(j)  is  the  total  energy  at  the  ,;'th  scale,  and  N(j)  is  the 
number  of  wavelet  coefficients  at  this  scale  (i.e.,  the  number  of  surface  coefficients). 
According  to  the  scale  designation  of  the  MATLAB  Wavelet  Toolbox, 

N(j)  =  \n(j  +  1).  (VI.7) 
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If  the  signal  x(t)  resides  within  the  scale  j  or  resides  within  the  scale  (j  + 1), 
with  the  same  total  energy,  we  have 

A(j  +  1)  =  2A(j) .  (VI.8) 

Thus,  there  will  be  a  3  dB  gain  in  the  energy  per  sample  per  each  increment  in  the 
scale  index. 

Table  1  summarizes  the  energy  distribution  obtained  with  wavelet  analysis 
using  Daubechies  wavelet  of  order  2  (Daub-2)  and  order  10  (Daub-10).  There  are 
three  input  vectors  with  64  samples  each,  the  first  vector  has  a  frequency  of  3/8 Fs, 
the  second  has  3/16F5,  and  the  third  has  3/32Fs,  where  Fs  is  the  sampling  frequency. 
This  means  that  the  first  input  vector  resides  within  the  first  scale,  the  second  vector 
resides  within  the  second  scale,  and  the  third  vector  resides  within  the  third  scale. 
We  conclude  the  following: 

•  The  total  energy  of  the  input  signals  is  distributed  among  the  scales.  The  sum 
of  the  total  energies  over  the  scales  is  slightly  less  than  the  input  signal  total 
energy  since  we  disregard  contribution  from  the  low  pass  section. 

•  The  proper  scale  (where  the  signal  resides)  has  the  greatest  share  of  the  total 
energy.  This  share  is  greater  (percentage  wise)  if  a  longer  wavelet  (e.g.,  Daub- 
10)  is  used  than  for  a  shorter  wavelet  (e.g.,  Daub-2). 

•  Energy  per  sample  at  the  proper  scale  is  larger  if  the  signal  resides  at  a  higher 
scale.  The  gain  factor  in  energy  per  sample  (if  the  signal  resides  within  scale  2 
rather  than  scale  1)  is  approximately  1.41  and  1.64  for  Daub-2  and  Daub-10, 
respectively.  We  note  that  ideally,  the  gain  factor  in  energy  per  sample  should 
be  2  per  scale  index. 
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Table  1:     Energy  Distribution  of  Daub-2  and  Daub- 10  Wavelets 


Measure 

Wavelet 

Scale 

/  =  (3/8)F, 

Signal  with 
/  =  (3/16)F. 

/  =  (3/32)Fs 

Max 
Coefficient 

Daub-2 

1 
2 
3 

1.3365 
0.5753 
0.3743 

0.6574 
1.6998 
0.3594 

0.4349 
1.0356 
1.8480 

Daub-10 

1 

2 
3 

1.2524 
0.3968 
0.1854 

0.3818 
2.0331 
0.3859 

0.3284 
0.5504 
2.8296 

Total 
Energy 

Daub-2 

1 
2 
3 

29.7687 
1.1313 
1.0495 

7.2565 

23.3282 

0.9492 

0.8486 

8.2775 

21.2601 

Daub-10 

1 

2 
3 

31.4508 
0.3980 

0.0875 

L6228 

29.9017 

0.3900 

0.4040 

1.9671 

29.0447 

Sample 

Energy 

(average 

per  sample) 

Daub-2 

1 
2 
3 

0.9201 
0.0628 
0.1049 

0.2199 
1.2960 
0.0949 

0.0257 
0.4599 
2.1260 

Daub-10 

1 

2 
3 

0.7671 
0.0133 
0.0036 

0.0396 
0.9967 
0.0163 

0.0099 
0.0656 
1.2102 
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2.       Scale  Identification  in  an  Ideal  Noise  Environment 

In  this  particular  simulation,  the  correlation  surface  is  approximated  by  si- 
nusoidal signals  embedded  in  ideal  white  Gaussian  noise.  The  maximum  coefficient, 
total  energy,  and  energy  per  sample  are  tested.  The  signal-to-noise  ratio  is  expressed 
as: 


SNR  = 


input  signal  power 


input  noise  variance 
Three  different  input  vectors  at  3Fs/8,  3FS/16,  and  3F5/32,  each  with  2048  samples, 

are  analyzed  using  Daub-2  and  Daub-10  wavelets.  SNR  values  are  10,  5,  2.5,  0,  —2.5, 

—5,  —10,  and  —15  dB.  Each  SNR  level  was  used  in  20  trials.  Results  are  shown  in 

Table  2  where  the  tabulated  SNR  is  the  smallest  SNR  at  which  each  measure  correctly 

designates  the  proper  scale  100%  of  the  time. 

Table  2:     Noisy  Observation 


Measure 

Wavelet 

Type 

Greatest 
Coefficient 

Total 
Energy 

Sample 
Energy 

Smallest  SNR  [dB] 

Daub-2 

2.5 

0 

-10 

Daub-10 

0 

-2.5 

-10 

These  results  show  that  the  energy  per  sample  achieves  the  lowest  SNR  value  com- 
pared to  the  other  two  measures.  Consequently,  this  measure  is  considered  the  most 
reliable  one  for  proper  scale  identification. 

3.       Hop-Scale  Pattern 

We  recall  that  scale  identification  means  assigning  a  scale  index  to  each  hop 
of  the  observed  signal.  The  hop  is  assigned  the  scale  whose  energy  per  sample  is  the 
greatest  among  the  values  of  the  other  scales;  hence,  it  is  called  the  proper  scale. 
Therefore,  a  sequence  of  hops  will  have  a  sequence  of  proper  scales,  also  called  the 
hop-scale  pattern. 
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Properly  selecting  a  heterodyne  and  a  sampling  frequency  will  result  in  proper 
choice  of  the  wavelet  scales.  The  number  of  scales  is  defined  by  the  number  of  octaves 
in  the  heterodyned  signal.  The  designation  of  the  scales  (scale  indices)  is  controlled 
by  the  selected  sampling  frequency,  as  each  scale  index  is  defined  in  terms  of  fractions 
of  the  sampling  frequency  as  shown  in  Figure  15.  As  an  example,  assume  we  have 
an  FH  signal  with  hopping  bandwidth  of  30-90  MHz.  If  the  signal  is  heterodyned  to 
10-70  MHz  and  sampled  at  280  MHz,  then  the  proper  wavelet  scales  are  S2,  S3,  and 
S4.  If  the  FH  signal  is  heterodyned  to  2-62  MHz  and  sampled  at  280  MHz,  then  the 
proper  wavelet  scales  are  S2,  S3,  S4,  S5,  and  S6. 

To  distinguish  the  FH  signal  from  other  modulation  types,  the  number  of 
scales  used  must  be  at  least  two.  The  hop-scale  pattern  will  then  be  comprised  of 
two  symbols  (i.e.,  two  scale  numbers).  For  other  modulations  the  hop-scale  pattern 
will  be  comprised  of  one  symbol  (i.e.,  one  scale  number).  For  the  MFSK  scheme  with 
a  carrier  frequency  in  the  above  mentioned  band,  the  bandwidth  of  the  modulated 
signal  is  typically  25  KHz  (i.e.,  the  channel  separation  for  the  wireless  sets).  Thus, 
the  MFSK  signal  bandwidth  can  not  span  an  octave  except  when  we  heterodyne  to 
a  very  low  frequency  (i.e.,  less  than  25  KHz).  Hop-scale  patterns  are  also  useful  for 
hop  frequency  extraction  from  wavelet  surfaces,  since  they  point  to  the  proper  scale 
where  most  of  the  signal  energy  is  concentrated.  Therefore,  frequency  extraction,  if 
done  on  that  wavelet  surface,  will  provide  correct  values. 

4.       Success  Rate 

Performance  of  scale  identification  is  evaluated  as  the  success  rate,  Pa.  We 
generate  known  hop-scale  patterns  and  obtain  the  percentage  of  the  correctly  identi- 
fied hop-scales.  Therefore,  Pid  is  defined  as: 

number  of  correct  hop-scales 
total  number  of  hops 

The  probability  of  error,  given  by  Pe  =  1  —  P^,  is  the  error  in  identifying  the  correct 
scale.  The  quality  of  scale  identification  depends  on  the  height  of  the  greatest  energy 
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per  sample  relative  to  the  other  sample  energies  from  other  scales.  Therefore,  a 
suggested  measure,  d^,  of  the  quality  of  scale  identification  is  defined  as: 

_  mM{A(j)})  -  mean({^(j)}) 

where  max(-)  is  the  maximum  value,  mean  (•)  is  the  mean  value,  and  var(-)  is  the 
variance.  The  term  dh  measures  the  distance  between  the  height  of  the  energy  per 
sample  at  the  proper  scale  to  its  average  values  over  all  scales,  and  is  expressed  in 
terms  of  the  standard  deviation  of  the  energy  per  sample. 

E.       EQUALIZATION  OF  THE  SPECTRAL  SHAPE  OF 

THE  INSTANTANEOUS  CORRELATION  FUNCTION 

Simulations  show  that  the  ICF  of  white  noise  (data  set)  has  a  triangular  type 
spectrum  in  the  direction  of  the  time  delay.  Figure  20  shows  theoretical  and  experi- 
mental frequency  spectra  taken  in  direction  of  the  time  delay  (i.e.,  transforming  the 
time  delay  to  frequency).  Thus,  the  energy  per  sample  for  all  wavelet  scales  need  to 
be  compensated  since  the  spectrum  is  colored. 

Figure  20  shows  that  the  slope  of  the  ICF  spectrum  is  1,  thus,  its  height,  q(j), 
at  the  middle  of  the  jth  scale  is: 

3(iV-l) 


q(j)  =  1  + 


2J+1 


where  N  is  the  number  of  ICF  data  points  in  the  direction  of  the  time  delay.  Note 
that  the  ICF  of  white  noise  is  created  by  multiplying  samples  of  the  white  noise  in 
the  time  domain.  Therefore,  the  resultant  FT  of  the  ICF  is  the  convolution  of  two 
FT  of  the  white  noise,  consequently,  we  obtain  the  triangular  spectral  (FT)  shape. 
As  a  result,  we  will  apply  an  equalizer,  whose  values  are  the  reciprocal  of  q{j),  to 
compensate  the  magnitude  distortion  at  all  scales. 

The  performance  of  the  scale  identification  is  evaluated  and  a  new  success  rate 
Pid  is  computed.  Results  are  presented  Chapter  VII. 


76 


F.       FREQUENCY  ESTIMATION 

In  Chapter  III,  we  related  the  WVD  of  the  signal  x(t)  to  the  wavelet  trans- 
form of  its  ICF.  The  relation  is  stated  in  Equation  III. 33  and  it  is  repeated  here  for 
convenience: 

/oo 
G*(sf)vx(tjy2*fldt, 
-00 

where  W*x(t;sJ)  is  the  wavelet  transform  of  the  ICF  of  x(t),  Vs(t,  f)  is  the  WVD 
of  x{t).  In  Chapter  2,  we  discussed  the  properties  of  WVD  and  concluded  that  the 
WVD  provides  the  means  to  estimate  the  signal  frequencies.  The  FT  of  the  wavelet 
surface  gives  a  bandpass  filtered  version  of  the  WVD  of  the  FH  signal.  Thus,  the  hop 
frequency  can  be  extracted  from  a  Fourier  transform  of  the  wavelet  surface,  in  the 
delay  direction,  over  the  main  diagonal  of  the  proper  diamond. 

1.       Frequency  Resolution 

The  Fourier  transform  of  the  wavelet  coefficients  can  be  used  for  spectral 
estimation  at  any  scale.  Consequently,  the  frequency  resolution  will  be  dependent 
on  the  Fourier  transformation.  The  frequency  resolution  (i.e.,  the  minimum  spacing 
between  two  resolved  narrow  band  components)  of  the  DFT  is  approximately  equal 
to  A/  =  Fs/N.  Thus,  at  any  given  scale  k,  the  number  of  data  points  N(k)  will  be 
related  to  the  number  of  input  data  points  N  by: 

N(k)  =  ffi. 

The  sampling  frequency  at  which  the  data  points  are  taken  is  scale  dependent,  i.e., 
Fs(k)  =  Fs/2k,  where  Fs  is  the  input  sampling  frequency. 

At  the  kth  scale,  both  the  number  of  data  points  (wavelet  coefficients)  and  the 
sampling  frequency  have  been  reduced  by  the  same  factor.  Consequently,  the  FT  of 
the  N(k)  data  points  has  a  frequency  resolution  given  by: 

*"*>-§§ -IM-  (vi.il) 

This  means  the  frequency  resolution  of  the  Fourier  transform  is  constant,  independent 
of  the  wavelet  scale. 
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2.  Success  Rate 

Performance  of  the  frequency  estimation  procedure  is  evaluated  in  terms  of  the 

success  rate,  Pf.  The  hop  frequency  is  considered  correctly  estimated  if  the  spectral 

peak  is  at  the  true  spectral  bin.  Pf  is  defined  as: 

number  of  correct  hop  frequencies 
total  number  of  hops 

The  probability  of  error  in  estimating  the  correct  hop  frequency  is  denoted  by  Pe, 

where  Pe  =  1  —  Pf.  The  hop  frequency  is  estimated  as  the  bin  number  corresponding 

to  the  peak  of  the  FT  over  a  specified  region  of  the  wavelet  surface.  Therefore,  the 

quality  of  the  frequency  estimation  depends  on  the  spectral  height  of  the  peak  relative 

to  the  average  background.  A  measure  of  the  quality  of  the  frequency  estimation  is 

given  by: 

d    =  max  (\X(f)\)-  mean  (\X(f)\) 

'   '  )/var(|X(/)|) 

where  X(f)  is  the  FT  over  the  wavelet  surface,  max  (•)  is  the  peak  value,  mean  (•)  is 

the  mean,  and  var  (•)  is  the  variance.  The  variable  df  measures  the  distance  between 

the  peak  value  of  X(f)  to  the  average  background  in  units  of  the  standard  deviation 

ofX(f). 

3.  Alternatives  for  Frequency  Estimation 

Given  the  correct  hop  start/stop  times,  the  hop  frequency  may  also  be  esti- 
mated directly  from  the  time  signal  or  from  the  ICF  surface.  To  extract  the  frequency 
from  the  original  time  signal  we  use  the  FFT  over  the  hop  length.  Recall  that  the 
FFT  is  a  matched  filter  for  sinusoidal  signals  in  white  Gaussian  noise.  Thus,  a  bet- 
ter performance  is  expected  in  comparison  to  the  nonlinear  processing  of  the  signal 
through  the  ICF  computation  and  the  wavelet  transformation.  The  performance  of 
frequency  estimation  from  the  signal  or  from  its  ICF  is  evaluated  using  the  already 
defined  measures  Pf  and  df.  We  note  that  Pfs  and  PfC  denote  the  success  rate  using 
the  original  time  signal  or  using  the  ICF  surface,  respectively.  Also,  dfs  and  dfC  de- 
note the  measure  of  the  frequency  estimation  quality  at  the  original  signal  or  at  the 
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ICF  surface,  respectively.  The  performance  of  PfS  and  PfC  are  evaluated  and  results 
will  be  given  in  Chapter  VII. 

G.      EXTRACTION  OF  HOP  TIMES 

Hop  timing  extraction  is  obtained  by  estimating  the  hop  start  and  stop  points. 
We  recall  that  ICF  surface  and  wavelet  surfaces  have  both  a  cellular  structure  consist- 
ing of  diamonds,  where  each  diamond  is  associated  with  a  specific  hop.  The  diamond 
edges  define  the  start/stop  point  of  the  hops.  The  diamond  width  corresponds  to 
the  hop  interval  7/j.  The  sides  (edges)  of  the  diamonds  of  hops  are  mutually  par- 
allel and  spaced  by  the  hop  interval  TV  Therefore,  we  obtain  the  distance  between 
two  sequential  intersections  of  the  diamonds  and  the  time  axis,  to  determine  the  hop 
start  (or  stop)  point.  This  tends  to  transform  the  transition  point  detection  into  an 
edge  detection  problem.  There  are  many  approaches  to  tackle  the  problem  of  hop 
timing  extraction;  in  the  following  section  we  address  one  technique  based  on  an  edge 
detection  operator. 

Compass  Operator 

Edge  detection  is  a  fundamental  problem  in  image  analysis  since  edges  help 
to  identify  objects.  Edges  are  characterized  by  an  abrupt  change  in  gray  level,  there- 
fore, edges  can  be  detected  using  the  derivative  (gradient)  operators  which  will  be 
maximum  in  the  direction  of  the  edges.  There  are  two  types  of  edge  operators,  the 
gradient  operators  and  the  compass  gradient  operators  (also  called  the  compass  oper- 
ators) [Ref.  62].  The  gradient  operator  measures  the  gradient  of  the  two-dimensional 
image  in  two  orthogonal  directions.  It  is  usually  applied  to  detect  edges  with  unknown 
directions.  The  compass  operator  measures  the  gradient  of  the  two-dimensional  im- 
age in  a  specific  direction.  It  is  used  to  detect  edges  with  pre-determined  directions. 
Since  the  wavelet  surface  has  edges  at  angles  of  §  and  *?  radians,  a  specific  com- 
pass operator  can  be  used  to  detect  the  edge  in  these  directions.  There  is  a  variety 
of  compass  operators.    They  differ  in  form,  depending  on  the  direction  of  the  edge 
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to  be  found.  An  example  of  the  compass  operator  is  the  NE  compass  (  NE  stands 
for  North-East).  It  has  a  3  x  3  matrix  NE  =  [0, 1, 1; -1,0, 1; -1, -1,0].  Another 
example  is  the  NW  compass  (NW  stands  for  North- West)  with  the  3x3  matrix 
NW  =  [1, 1, 0;  1, 0,  —1, 0,  —1,  —1].  The  compass  operator  is  applied  to  the  upper  half 
of  the  wavelet  surfaces  to  potentially  detect  the  diamond  edges.  An  array  of  the  NE 
compass  operator  is  shown  in  Figure  21.  The  compass  array  is  used  to  scan  the  sur- 
face from  left  to  right.  All  of  the  contributions  are  summed  according  to  the  weights 
of  the  3x3  matrix.  The  maximum  value  is  extracted  to  define  the  point  where  the 
compass  array  matches  an  edge.  To  make  our  data  applicable  to  compass  operations 
we  need  to  take  care  of  the  negative  portions  of  the  surfaces,  which  is  done  by  adding 
in  a  positive  number  equals  in  magnitude  to  the  smallest  (i.e.,  the  most  negative) 
surface  value.  Results  are  presented  in  Chapter  VII. 
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Figure  14.  Functional  description  of  the  processing  scheme. 
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Figure  15.  MATLAB  wavelet  scale  designation  and  computation  tree. 
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Figure  16.  Wavelet  surfaces  of  FH  signals  using  Daub-3  at  10  dB  (scale  1,  2,  and  3). 
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Figure  17.  Wavelet  surfaces  of  FH  signals  using  Daub-3  at  10  dB  (scale  4  and  5). 
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Figure  18.  Wavelet  surfaces  of  FH  signals  using  Daub-3  at  3  dB  (scale  1,2, and  3), 
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Figure  19.  Wavelet  surfaces  of  FH  signals  using  Daub-3  at  3  dB  (scale  4  and  5). 
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Figure  20.  Theoretical  and  simulated  spectrum  of  the  ICF  of  white  noise. 
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Figure  21.  The  compass  line  over  the  upper  half  of  the  wavelet  surface. 
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VII.         SIMULATIONS  AND  RESULTS 

A.  INTRODUCTION 

In  this  chapter,  we  will  provide  experimental  results  of  the  processing  scheme 
introduced  in  Chapter  VI.  Section  VII. B  shows  the  results  obtained  by  visual  inspec- 
tion of  wavelet  surfaces  as  they  pertain  to  the  identification  of  FH  signals.  Wavelet 
surfaces  from  other  digital  modulation  schemes  are  compared  with  those  resulting 
from  the  FH  signals.  Section  VII.C  provides  results  of  hop-scale  identification.  Sec- 
tion VII.D  provides  results  of  hop  frequency  extraction  and  compares  results  obtained 
using  frequency  estimation  from  wavelet  surfaces  and  from  the  original  time  signal. 
Section  VILE  provides  the  results  of  hop  start/stop  times  estimation  obtained  by 
using  the  compass  operator  over  the  wavelet  surface. 

B.  VISUAL  INSPECTION 

We  carried  out  an  opinion  test  to  identify  the  FH  scheme  by  examining  the 
wavelet  surfaces.  Ten  participants  were  involved,  each  one  was  asked  to  identify  the 
diamond  structure  of  the  FH  signal  over  the  wavelet  surfaces  at  all  pertinent  scales  and 
for  all  hops.  Two  types  of  wavelet  were  used;  the  Morlet  wavelet  and  the  Daubechies 
wavelet  of  order  3.  Both  wavelets  were  used  in  their  real  and  complex  form.  Four 
SNR  values  were  used;  10,  6,  3  and  0  dB.  Four  different  surface  representations 
were  used;  the  real  part,  imaginary  part,  magnitude,  and  phase.  To  avoid  biasing 
and  preconception  all  participants  started  to  identify  the  surfaces  resulting  from  the 
lowest  SNR  value  and  then  the  higher  SNR  values  consequently.  Detailed  scoring 
tables  and  the  scoring  code  are  given  in  Appendix  B.  Only  scoring  tables  belonging 
to  10  and  3  dB  are  given  in  this  section  because  of  their  significance  to  the  final 
conclusion  of  the  opinion  test.  10  dB  is  the  highest  tested  SNR  value  while  3  dB 
is  the  minimum  SNR  which  provided  the  minimum  acceptable  identification  score  of 
2  (according  to  the  scoring  code).   In  Appendix  B,  each  scoring  table  is  dedicated 
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for  a  specific  wavelet  type  at  a  specific  SNR  value  and  contains  the  scores  given  for 
different  representations  at  different  wavelet  scales.  "Real"  stands  for  representing 
the  wavelet  surface  by  the  contour  plot  of  its  real  part,  while  "Imag."  ,  "Abs."  and 
"Angle"  stand  for  representing  the  surface  using  the  contour  plots  of  imaginary  part, 
magnitude  and  phase,  respectively.  The  fractions  in  numbers  listed  in  the  tables  came 
from  averaging  the  ten  scores  given  by  the  ten  participants. 

Table  3  shows  the  summary  of  scores  obtained  from  different  wavelet  types  at 
SNR  values  of  10  and  3  dB  when  representing  the  wavelet  surfaces  by  their  real  part. 
Eventually,  the  real  part  and  the  imaginary  part  outperformed  other  representations. 
The  ratings  in  Table  3  ranges  from  0.2  to  1,  where  1  represents  perfect  identification 
of  the  hop  diamonds  at  their  proper  positions  and  0.2  represents  just  distinction  of 
the  hop  diamonds  from  the  background  noise.  Test  signals  were  designed  to  reside 
within  the  five  wavelet  scales  given  in  the  table.  Comparing  different  scoring  tables 
in  Appendix  B,  it  becomes  apparent  that  the  real  part  or  imaginary  part  provided 
the  best  surface  representation  for  the  purpose  of  visual  inspection.  The  magnitude 
or  phase  provided  a  very  poor  representation.  We  noted  that  real  valued  wavelets  did 
a  slightly  better  job  than  the  complex  valued  wavelets  (at  least  for  the  wavelet  types 
used  in  the  simulations). 

Thus,  this  visual  opinion  test  indicates  that: 

(1)  The  FH  signal  can  be  identified  over  the  wavelet  surfaces  by  its  cellular  struc- 
ture which  is  dominant  at  the  proper  scale.  That  is,  each  scale  will  emphasize 
the  hops  which  belong  to  the  scale  and  attenuate  other  hops.  The  (diamond) 
cellular  structure  provides  an  estimate  for  the  hop  start/stop  times. 

(2)  The  FH  signal  can  be  well  identified  at  3  dB  SNR  and  above. 

(3)  The  real  part  or  imaginary  part  provides  the  best  surface  representation  for 
the  purpose  of  visual  inspection. 

(4)  The  real  form  of  the  wavelet  function  provides  better  surface  representation 
than  the  complex  form  for  the  wavelets  used  in  simulations.  However,  other 
types  of  wavelets  may  perform  differently. 
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(5)  The  Morlet  wavelet  has  better  performance  than  the  Daubechies  wavelet  es- 
pecially at  higher  scales.  This  is  because  our  Morlet  wavelet  has  a  narrower 
bandwidth  than  that  obtained  with  Daubechies  wavelet,  and  the  Daubechies 
wavelet  results  in  a  decimated  transform  which  reduces  the  number  of  coeffi- 
cients at  higher  scales. 

(6)  Other  modulation  schemes  such  as  ASK,  PSK,  and  MFSK  will  have  wavelet 
surfaces  residing  only  at  one  scale.  Plots  of  these  wavelet  surfaces  are  given  in 
Appendix  C. 
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Table  3:  Summary  of  identification  scores  for  real  part  representation. 


Wavelet 
Type 

SNR 
[dB] 

Scales 

SI 

S2 

S3 

S4 

S5 

Real 

10 

1.0 

1.0 

1.0 

0.95 

0.95 

Morlet 

3 

1.0 

1.0 

1.0 

0.95 

0.9 

Complex 

10 

1.0 

1.0 

1.0 

0.9 

0.75 

Morlet 

3 

0.7 

0.8 

0.9 

0.5 

0.4 

Real 

10 

1.0 

1.0 

0.8 

0.2 

0.2 

Daub-3 

3 

0.8 

0.7 

0.5 

0.2 

0.2 

Complex 

10 

1.0 

1.0 

1.0 

0.3 

0.2 

Daub-3 

3 

0.6 

0.5 

0.7 

0.4 

0.2 

Scores  vary  between  0.2  and  1,  where 
1:  perfect  identification  of  hop  diamonds. 


•  0.2:  just  distinction  of  hops  from  background  noise. 
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C.      SCALE  IDENTIFICATION 

The  proposed  processing  scheme  is  used  to  extract  hop  start/stop  times,  the 
hop-scale  pattern,  and  the  hop  frequency.  First,  we  investigate  scale  identification 
and  hop  frequency  extraction  assuming  correct  hop  timing  is  available.  Extraction 
of  hop  start/stop  times  will  be  considered  later.  For  scale  identification  performance 
we  carried  out  simulations  using  the  following  data: 

(1)  Signal  pattern  length:  5  hops. 

(2)  Wavelet  scales:  SI,  S2,  S3,  and  S4. 

(3)  Wavelet  types:  Daub-2,  4,  and  8. 

(4)  SNR:  from  10  to  -10  dB. 

(5)  Number  of  experiments:  250  per  wavelet  per  SNR  value. 

The  hop  frequencies  of  FH  signals  are  selected  to  generate  the  hops  according 
to  a  known  scale  pattern  taking  in  consideration  the  doubling  effect  of  frequencies  by 
computing  the  ICF  surface.  The  test  pattern  is  selected  of  four  scales;  SI,  S2,  S3  and 
S4.  The  wavelet  surfaces  are  generated  from  the  ICF  surfaces  at  those  relevant  scales. 
The  total  energy  of  each  hop  at  each  scales  are  computed  and  the  energy  per  samples 
are  obtained  by  dividing  the  total  energy  by  the  number  of  wavelet  coefficients  at 
each  scale.  Over  each  hop,  the  scale  with  the  greatest  energy  per  sample  is  estimated 
as  the  proper  scale.  The  resultant  estimated  hop  pattern  is  compared  to  the  known 
hop  pattern  and  the  probability  of  correct  identification  is  computed.  We  compute 
the  energy  contained  in  each  hop  over  the  diamond  corresponding  to  that  hop.  We 
disregard  energy  contained  outside  the  diamond  since  it  is  contributed  by  different 
hops.  To  avoid  bias  from  the  colored  noise  of  the  ICF  surfaces  an  equalization  must 
be  performed  to  the  resultant  energy  per  samples  at  all  pertinent  scales  before  using 
them  in  deciding  upon  the  proper  hop  scale.  Although,  the  equalization  process  is 
based  on  the  theoretical  spectral  shape  given  in  Figure  20,  the  resultant  equalized 
energy  per  samples,  for  ideal  input  white  Gaussian  noise,  do  not  have  the  same  value 
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at  all  scales.  Table  4  presents  averages  and  variances  obtained  for  the  equalized 
energy  per  samples  for  Daubechies  of  order  2,  4,  and  8.  We  note  that  the  average 
values  for  Daub-2  and  Daub-8  are  monotonically  increasing  with  the  scale  index.  For 
Daub-4  they  are  monotonically  increasing  for  the  scales  S1,S2  S3  and  S5  while  they 
have  a  dip  at  scale  S4.  These  different  trends  suggest  that  different  wavelets  respond 
differently  to  the  theoretical  equalization  of  the  spectral  shape  of  the  ICF.  Therefore 
the  correction  weights  of  the  energy  per  sample  depend  on  the  wavelet  used.  In  these 
simulations,  results  are  given  for  the  proper  scale  identification  after  equalization  and 
correction. 

Figures  22  to  24  show  the  success  rate,  Pid  obtained  for  Daubechies  wavelets 
of  order  2,  4,  and  8  at  scales  SI,  S2,  S3,  and  S4.  For  the  performance  of  scale 
identification  we  consider  the  minimum  SNR  at  which  P^  is  1  as  the  figure  of  merit. 
Over  all  tested  scales  the  success  rate,  P^  assumes  the  value  of  1  at  different  minimum 
input  SNR  values.  This  is  a  function  of  the  order  of  the  wavelet  and  the  scale.  Figure 
22  shows  that  the  performance  of  Pa  obtained  from  Daub-2  has  achieved  the  value 
of  1  at  SNR  equals  —1  dB  at  most  of  the  scales,  hence  ,  —1  dB  is  considered  the 
minimum  SNR  value  for  Daub-2.  The  minimum  SNR  value  for  Daub-4  and  Daub-8 
is  —2  dB  at  most  of  the  scales  as  shown  in  Figures  23  and  24.  We  note  also  that 
for  Daub-8,  Pa  assumes  the  value  of  0.9  for  most  of  the  scales  at  an  SNR  of  —3  dB. 
Results  show  that: 

(1)  For  long  wavelets  (Daub-8)  the  probability  of  correct  scale  identification  is  1.0 
at  SNR  equals  to  -2  dB,  while  it  is  0.9  at  SNR  equals  -3  dB. 

(2)  Longer  wavelets  perform  better  than  shorter  ones. 

(3)  The  performance  at  higher  scales  is  mostly  better  than  the  performance  at 
lower  scales  in  terms  of  minimum  SNR  at  which  Pa  is  1.  The  exception  of 
that  case  is  the  performance  at  the  scale  SI  which  can  be  justified  by  the 
imperfect  equalization  of  the  ICF  spectral  shape. 


94 


Table  4:  Mean  and  standard  deviation  of  energy  per  sample. 


Mean 

SI 

S2 

S3 

S4 

S5 

Daub-2 

1.5688 

1.7557 

1.9479 

2.3165 

2.6086 

Daub-4 

1.6303 

1.8432 

2.0218 

2.5585 

2.1149 

Daub-8 

1.7173 

1.9979 

2.2725 

2.4776 

2.7525 

Standard 
Deviation 

SI 

S2 

S3 

S4 

S5 

Daub-2 

0.3419 

0.3763 

0.6503 

0.9542 

1.3176 

Daub-4 

0.3507 

0.4103 

0.6891 

1.1328 

1.2776 

Daub-8 

0.3665 

0.4520 

0.7874 

1.1063 

1.5845 

Theoretically,  wavelet  surfaces  at  higher  scales  have  higher  SNR  values  than 
those  at  lower  scales  due  to  the  higher  processing  gain  at  higher  scales.  Consequently, 
the  decay  in  the  Pjd  performance,  at  low  input  SNR,  is  slower  at  higher  scales  than 
the  decay  at  lower  scales.  But,  we  noted  that  the  decay  in  Pa  is  not  consistent 
at  different  scales,  especially  at  scale  SI,  which  shows,  mostly,  the  slowest  decay 
compared  to  other  scales.  This  anomaly  may  come  from  many  interacting  effects  due 
to  the  imperfect  equalization  and  the  non-ideal  wavelet  filters  at  different  scales. 

D.      FREQUENCY  EXTRACTION 

We  carried  out  simulations  to  evaluate  the  performance  of  frequency  estimation 
using  the  data  from  Section  VII.C. 
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The  hop  frequencies  are  estimated  by  taking  the  FT  of  the  vector  of  wavelet 
coefficients  located  at  the  center  of  the  diamond  in  the  direction  of  the  time  delay. 
The  bin  corresponding  to  the  peak  value  represents  the  estimated  hop  frequency.  The 
estimated  hop  frequency  is  then  compared  to  the  true  hop  frequency  and  the  proba- 
bility of  correct  frequency  extraction  (the  success  rate)  is  computed.  The  estimated 
frequency  is  considered  correct  if  the  estimation  error  is  less,  in  percentage  of  the  true 
frequency,  than  -h,  where  TV  is  the  length  of  the  vector  of  the  wavelet  coefficients. 

Figures  25  to  27  plot  the  success  rate  Pf  obtained  for  different  wavelets  at 
different  scales.  Figures  28  to  30  display  the  corresponding  values  of  the  distance  df. 

We  consider  the  minimum  SNR,  at  which  Pf  is  1,  as  the  figure  of  merit  for 
the  performance  of  frequency  extraction.  Over  all  tested  scales  the  success  rate  of 
frequency  estimation,  Pf,  assumes  the  value  of  1  at  different  minimum  input  SNR 
values. 

Figure  25  shows  that  the  Pf  value  for  Daub-2  is  1  at  SNR  equals  0  dB  at  most 
of  the  scales,  hence,  0  dB  is  considered  the  minimum  SNR  value  for  Daub-2.  The 
minimum  SNR  value  for  Daub-4  and  Daub-8  is  also  0  dB  at  most  of  the  scales,  as 
shown  in  Figures  26  and  27. 

For  Daub-2,  Daub-4  and  Daub-8  the  minimum  input  SNR  is  0  dB  at  scales 
SI,  S2,  and  S3.  At  a  success  rate  Pf  of  0.9,  the  minimum  SNR  equals  —2  dB  over 
most  scales.  The  distance  df  is  about  4  times  the  standard  deviation  of  the  wavelet 
surface  data  for  all  tested  scales  as  shown  in  Figures  28  to  30.  Values  of  the  success 
rate  Pf  and  the  distance  df  shows  that  the  hop  frequency  can  be  reliably  extracted 
from  the  wavelet  surfaces  at  an  SNR  equals  0  dB  or  better  using  only  one  FT  at  the 
center  of  the  hop  diamond  area  in  the  direction  of  the  time  delay. 
Alternatives  for  frequency  estimation 

Hopping  frequencies  may  also  be  estimated  from  the  original  signal  or  from  the  ICF. 
Figure  31  plots  the  performance  Pfs  and  Pfc  obtained  from  the  time  signal  and  from 
the  ICF,  respectively,  at  different  SNR  values,  assuming  exact  estimates  of  hopping 
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start/stop  times  are  available.  Figure  32  plots  the  distance  dfS  and  djc  obtained 
from  the  time  signal  and  from  the  ICF,  respectively,  at  different  SNR  values.  The 
frequency  estimation  success  rate  from  the  time  signal,  PfS,  is  1  at  SNR  value  of 
—5  dB.  The  frequency  estimation  success  rate  from  the  ICF,  PfC,  is  1  at  SNR  value 
of  0  dB.  The  distances  dfS  and  dfc,  obtained  from  the  original  signal  and  from  the 
ICF  are  6  and  7,  respectively.  This  shows  that,  assuming  exact  estimates  of  hopping 
start/stop  times  are  available,  hop  frequencies  may  be  estimated  by  processing  the 
original  signal  at  lower  SNR  values  than  can  be  achieved  using  the  wavelet  surfaces. 
The  benefit  obtained  by  analyzing  the  ICF  surface  by  wavelet  analysis  is  significant 
in  case  of  unknown  hop  start/stop  times.  Also,  less  computations  are  needed,  in  case 
of  estimating  the  frequency  from  wavelet  surfaces,  due  to  applying  one  FT  per  hop, 
using  few  coefficients. 

E.       EXTRACTION  OF  HOP  TIMES 

This  section  presents  hopping  time  extraction  results  obtained  using  the  com- 
pass operator  discussed  earlier  in  section  VI. G.  The  wavelet  surface  is  represented 
by  its  upper  half  plane  resulting  in  a  triangular  shape  instead  of  the  familiar  hop 
diamond  shape.  The  line  compass  operator  is  used  over  the  surface  from  left  to  right 
and  the  peak  value  of  the  resultant  contribution  is  located  over  the  time  axis.  The 
location  of  the  peak  value  points  to  the  estimated  hop  starting  time.  The  difference 
between  the  true  and  the  estimated  starting  time  is  considered  the  estimation  error 
and  is  evaluated  in  terms  of  points  over  the  time  axis.  Each  SNR  value  is  used  in  20 
trials.  Figures  33-35  show  the  results  of  the  estimation  error  (Er)  at  different  SNR 
values  for  different  Daubechies  wavelet  lengthes.  The  mean-square  error  (MSE)  and 
standard  deviation  (SD)  of  the  estimation  error  are  also  given.  Results  show  that 
the  hop  start/stop  times  can  be  extracted  within  accuracy  of  ±7  points  (out  of  the 
128-point  hop  interval)  at  SNR  values  of  6  dB  or  better. 
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Figure  22.  Pid  for  Daub-2. 
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Figure  23.  Pid  for  Daub-4. 
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Figure  24.  Pid  for  Daub-8. 
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Figure  25.  Pf  for  Daub-2. 
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Figure  26.  Pf  for  Daub-4. 
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Figure  27.  Pf  for  Daub-8. 
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Figure  28.  The  spectral  distance  dt  for  Daub- 2. 
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Figure  29.  The  spectral  distance  dj  for  Daub-4. 
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Figure  30.  The  spectral  distance  df  for  Daub-8. 
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Figure  31.  Pjs  using  the  signal  and  PfC  using  the  ICF. 
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Figure  32.  dfs  using  the  signal  and  dfc  using  the  ICF. 
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Figure  33.  Hop  timing  estimation  error  for  Daub-2. 
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Figure  34.  Hop  timing  estimation  error  for  Daub-4, 
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Figure  35.  Hop  timing  estimation  error  for  Daub-8. 
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VIII.         CONCLUSIONS  AND 
RECOMMENDATIONS  FOR  FUTURE  WORK 

A.      CONCLUSIONS 

Wavelet  analysis  of  correlation  functions  is  a  new  area  which  can  be  used 
to  intercept  communication  signals  embedded  in  noise.  This  work  aimed  at  apply- 
ing wavelet  analysis  to  the  instantaneous  correlation  function  to  identify  frequency 
hopped  signals. 

The  pure  FH  is  the  most  popular  scheme  for  frequency  hopped  signals.  The 
instantaneous  correlation  function  (ICF)  of  the  complex-valued  pure  FH  signal  was 
shown  to  have  a  cellular  (diamond)  structure,  where  each  hop  contributes  to  one  main 
diamond.  Inside  this  diamond,  the  ICF  has  a  single  complex  exponential  component 
representing  the  hop  frequency  in  the  delay  direction.  We  showed  that  the  diamond 
intersections  with  the  time  axis  are  the  hop  start/stop  times  while  the  width  of  the 
diamond  is  the  hop  interval. 

The  wavelet  transform  of  the  ICF  surface  results  into  a  number  of  surfaces 
each  at  a  specific  wavelet  scale  (called  the  wavelet  surface).  The  wavelet  surface  at 
any  scale  emphasizes  the  hops  which  reside  in  it  and  attenuates  other  hops. 

We  addressed  the  interception  problem  using  two  possible  approaches.  In  the 
first  approach,  we  visually  inspected  the  wavelet  surfaces  to  identify  and  classify  the 
structure  of  the  FH  signal  and  to  obtain  a  rough  estimate  for  the  hop  time  interval. 
In  the  second  approach,  we  applied  the  processing  scheme  which  can  also  be  used  to 
automate  the  interception  task.  The  proposed  processing  scheme  estimates  the  hop 
start/stop  times,  the  hop-scale  pattern,  and  the  hop  frequency.  The  extraction  of  the 
hop  start/stop  times  was  addressed  using  an  edge  detection  approach,  by  applying 
a  compass  operator,  which  is  a  well  known  technique  in  the  image  processing  area. 
The  hop-scale  pattern  is  obtained  by  applying  an  energy  analysis. 

The  frequency  of  each  hop  can  be  extracted  from  the  wavelet  surface  at  the 
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proper  scale,  or  as  an  alternative,  from  the  original  time  data  using  the  time  param- 
eters. 

Visual  inspection  of  the  wavelet  surfaces  showed  that  the  FH  signal  can  be 
identified  from  the  wavelet  surfaces  for  an  input  SNR  of  3  dB  and  above.  Other 
modulation  schemes  such  as  ASK,  PSK,  and  MFSK  will  have  significant  wavelet 
surfaces  at  one  scale  only  since  their  frequency  bandwidth  does  not  span  more  than 
one  wavelet  scale.  Hop  timing  extraction  showed  that  the  hop  start/stop  times  can 
be  extracted  within  accuracy  of  ±5%  for  short  wavelets  (Daub-2)  at  SNR  of  6  dB  or 
better. 

For  the  hop-scale  identification,  we  adopted  the  energy  per  sample  as  a  measure 
for  the  proper  scale,  assuming  exact  estimates  of  hop  start  /stop  times  are  available. 
Results  showed  that  the  probability  of  correct  scale  identification  is  1.0  at  an  input 
SNR  of  -2  dB  for  long  wavelets  (Daub-8),  and  it  is  0.90  at  an  input  SNR  of  -3 
dB.  The  performance  of  longer  wavelets  is  better  than  that  of  shorter  ones  since 
longer  wavelets  have  better  spectral  energy  concentration  than  shorter  ones.  For 
hop  frequency  extraction,  the  success  rate  of  frequency  estimation  from  the  wavelet 
surfaces  showed  that  the  probability  of  correct  frequency  extraction  is  1.0  at  an  input 
SNR  of  0  dB  and  above. 

Results  showed  that  the  minimum  SNR  required  for  automating  the  intercep- 
tion task  is  6  dB.  However,  a  better  performance  of  frequency  estimation  is  obtained 
at  lower  SNR  value  (—5  dB)  by  processing  the  original  time  signal,  using  the  hop 
time  information. 

B.      RECOMMENDATIONS  FOR  FUTURE  WORK 

For  future  extension  of  this  work  we  recommend  the  following: 

(1)  Addressing  the  automatic  recognition  of  the  cellular  structure  of  the  FH  signal 
over  the  wavelet  surfaces.  There  are  two  topics  in  the  image  processing  area 
which  may  be  helpful  in  this  task;  the  automatic  recognition  of  regions  of 
interest,  and  automatic  target  recognition  [Ref.  63]. 
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(2)  Improving  the  performance  of  the  hop-scale  identification  at  lower  SNR  levels. 
In  addition,  it  is  suggested  to  readdress  the  equalization  of  the  spectrum  of 
the  ICF  surfaces. 

(3)  Investigating  other  wavelet  types,  and  the  use  of  other  definitions  for  the  in- 
stantaneous correlation  function. 

(4)  Combining  information  from  different  wavelet  surfaces  to  improve  parameter 
estimation. 
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APPENDIX  A.  THE  INSTANTANEOUS 
CORRELATION  FUNCTION  (ANALYTIC 

APPROACH) 


Let  the  complex  FH  signal  be  defined  as: 


N 


s(t)  =  exp    jJ2{unt  +  0n) 


71=1 


(A.l) 
te{n}, 


where  un  is  the  hopping  frequency  (in  radians)  and  9n  is  the  carrier  phase  chip  of  the 
nth  hop.  N  is  the  total  number  of  observed  hops.  The  shorthand  notation  t  €  {n} 
stands  for  the  condition: 

t0  +  (n  -  \)Th  <t<tQ  +  nTh,  (A.2) 

i.e.,  it  confines  the  value  of  the  time  t  corresponding  to  the  nth  hop  duration.  Th 
denotes  the  hop  interval  and  t0  denotes  an  initial  misalignment.  The  resulting  in- 
stantaneous correlation  function (ICF)  is  defined  in  terms  of  the  signal  values  at  time 
(t  +  r/2)  and  (t-r/2). 

If  "£"  is  restricted  over  the  nth  hop,  the  sum  term  (t  ±  r/2)  may  stay  over  the 
nth  hop  or  may  extend  to  (n  ±  1)  hop.  Consider  a  positive  "r",  i.e.,  0  <  r  <  Th, 
and  let  t  be  restricted  to  the  first  half  of  the  nth  hop  and  let  to  =  0,  then,  t  will  be 
restricted  to 

{n-l)Th<t<(n-^jTh,  (A.3) 

and  this  is  denoted  as: 

te{\n}'  (A-4) 

and  implies  that  if  t  is  restricted  to  the  1st  half  of  the  nth  hop  s(t  +  r/2)  will  remain 
within  the  nth  hop  since  0  <  r  <  TV  Therefore, 

s(t  +  ^)=exp(jJ2i^n(t  +  ^)  )).  (A.5) 

V  2J  \    «=1  I  V  2J   *+(T/2)6{n},  «6{(l/2)n}j; 
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If  (t  -  r/2)  stays  within  the  nth  hop  or  the  (n  -  l)st  hop,  then, 


st 


=    exp  (  Y^  I  un  1 t  - 


t-(r/2)G{n>,  t€{(l/2)n} 


f~  2 


(A.6) 
t-(r/2)€{n-l},  t€{(l/2)n}  J  / 

The  condition  t  €  {l/2n}  is  common  for  Equations  A. 5  and  A.6.  Therefore  we  take 
it  from  the  equation  body  and  place  it  aside  as  follows: 


R%r)    =    expjf£  Ln  (t  +  ^ 


t+(r/2)e{n}  V  2 


T 


t-(r/2)€{n} 


-CJn_i   ( t  -  - 


te 


IH 


(A.7) 

t-(r/2)€{n-i} . 

The  resultant  terms  of  this  summation  depend  on  the  validity  of  the  indica- 
tor functions,  i.e.,  terms  to  be  combined  together  must  have  common  (intersection) 
regions  of  their  indicator  functions.  Thus, 


R[{t,r)    =    expjljrlun(t  +  ^-un(t-Q 


t+(r/2)€{n},  *-(t/2)g{»} 


t+(r/2)€{n},  t-(r/2)€{n-l> 


(A.8) 


Let  us  examine  the  regions  satisfying  the  following  conditions: 

T  T 

Condition  1  :       t  +  —  €  {n},  t  -  -  €  {n}, 


Condition  2 


*  +  §€{n},  *-^e{n-l}, 


for     t  e  I  —n 


(A.9) 


,r  >0. 

The  region  where  Condition  1  is  valid  can  be  defined  from  Equations  A. 2  and 
A.  10  for  (t0  =  0)  as: 

(n-l)Th<t  +  ^<Th    and    (n-  l)Th  <  t-  ^  <  nTh,  (A.10) 
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where 

(n-l)Th<t<  (n-^\Th    and    Th  >  r  >  0.  (A.ll) 

Equation  A. 10  and  A.ll  can  be  plotted  according  to  the  following  set  of  inequalities 
(take,  e.g.,  n  =  1): 

*  +  5<Tfcl  *+,|>0, 
*-|<T&,  t-|>0, 


t  >  0,  t  <  (|)  Th, 


0<  r<Th 


(A.12) 


An  illustration  of  the  corresponding  set  of  equalities  is  shown  in  Figure  36. 
Therefore,  condition  1  is  valid  over  the  triangle  ABC.  Similarly,  condition  2  is 
defined  by  the  following  set  of  inequalities: 

t  +  ^<Th,  t  +  ^>0,  t-^>-Th,  t-^<0,  (A.13) 

and  0  <  t(Th/2),  0  <  r  <  T^.  Condition  2  is  valid  over  the  triangle  ADC,  therefore 
R[{t,r)    =    expjfx>n  (*  +  !)-"»(*-!) 

+wn  f  t  +  -  J  -  wn_i  f  t  -  -J 


Condition  1 


Condition  2 


or 


N 


R\t,r)  =  expj     ^u)nr 


.n=l 


+ 


Condition  1 


((Jn  -  CJn-i)t  +  ((Jn  +  wn_i)- 


,)■ 


Condition  2y 

(A.14) 


where  0  <  r,  t  e  {(l/2)n}. 

Note  for  all  hops,  condition  1  and  condition  2  valid  areas  will  be  as  shown  in  Fig.  37 
where  the  areas  of  valid  condition  1  are  denoted  by  1  and  the  areas  of  valid  condition 
2  are  denoted  by  2. 
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For  the  case  of  negativ  r,  while  t  is  still  in  the  1st  half  of  the  nth  hop,  we  make 
use  of  the  conjugate  symmetry  property  of  the  ICF,  this  will  conclude  Condition  3 
and  Condition  4  as  shown  in  Figure  38.  When  t  is  within  the  second  half  of  the  nth 
hop,  similar  procedure  will  lead  to  the  other  4  conditions  (  i.e.,  Condition  5,  6,  7,  8). 
These  areas  are  shown  in  Figure  39. 

Let  us  now  summarize  our  observations  from  the  derived  formulai  of  the  ICF 
by  introducing  the  doubly  indexed  function  Rm,n(t,T).  Rm,n(t,  r)  will  represent  the 
ICF  of  the  FH  signal  over  the  t  —  r  plane  inside  and  outside  the  diamond  of  each  hop. 
Rm,n{t,T)  will  give  the  ICF  expression  once  we  assign  its  indices  (m,  n).  The  indices 
are  shown  in  Figure  9.  Rm,n(t,T)  is  given  by 

#m,n(*,  r)  =  expj  J  (um  -  un)t  +  (um  +  u)n)-  J  ,  (A.15) 

where  m  and  n  are  the  indices  of  the  two  hops. 
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>►  t 


Figure  36.  Geometric  definition  of  Condition  1  and  Condition  2. 
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Figure  37.  Areas  of  Condition  1  and  Condition  2. 
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Figure  38.  Areas  of  Condition  3  and  Condition  4. 
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Figure  39.  Areas  of  Condition  1,2,. ..,8. 
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APPENDIX  B.  SCORES  OF  OPINION  TEST 

This  appendix  includes  the  tables  of  scores  of  the  opinion  test  carried  out  for 
visual  inspection  and  identification  of  the  FH  signal.  Only  one  table  is  shown  for  each 
wavelet  type.  Each  table  is  the  average  of  10  corresponding  tables  originated  from  10 
experiments.  This  justifies  the  fractions  in  the  tabulated  scores. 

Table  B.l:  Wavelet:  Complex  Daub-3,  SNR:10  dB 


s 

Real 

Imag. 

Angle 

Abs. 

1 

5 

5 

4.25 

1.25 

2 

5 

5 

4.75 

1 

3 

4 

5 

2 

1 

4 

1 

1.5 

1 

1 

5 

1 

1 

1 

1.25 

Table  B.2:  Wavelet:  Real  Daub-3,  SNR:10  dB 


s 

Real 

Imag. 

Angle 

Abs. 

1 

5 

5 

2.75 

1 

2 

5 

5 

5 

1 

3 

5 

5 

2 

1.25 

4 

1.5 

1 

1 

2.5 

5 

1 

1 

1 

1.5 

5  —  sharply  obvious  with  sharp  boundaries. 

4  —  obvious. 

3  —  fair. 

2  —  distinguishable  from  background  noise. 

1  —  not  distinguishable  from  background  noise. 
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Table  B.3:  Wavelet:  Complex  Morlet,  SNR:10  dB 


s 

Real 

Imag. 

Angle 

Abs. 

1 

5 

5 

4.5 

1 

2 

5 

5 

4.5 

1 

3 

5 

5 

4.5 

1 

4 

4.25 

4.75 

1 

3.5 

5 

3.75 

3.75 

1 

2.75 

Table  B.4:  Wavelet:  Real  Morlet,  SNR:10  dB 


s 

Real 

Imag. 

Angle 

Abs. 

1 

5 

5 

4.5 

4.5 

2 

5 

5 

4.5 

4.5 

3 

5 

5 

4.5 

4.75 

4 

4.75 

4.75 

1 

3.5 

5 

4.75 

3.75 

1 

3.75 

5  —  sharply  obvious  with  sharp  boundaries. 

4  —  obvious. 

3  —  fair. 

2  —  distinguishable  from  background  noise. 

1  —  not  distinguishable  from  background  noise. 
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Table  B.5:  Wavelet:  Complex  Morlet,  SNR:3  dB 


s 

Real 

Imag. 

Angle 

Abs. 

1 

3.5 

3.25 

3.5 

4 

2 

4 

3.5 

4 

3.5 

3 

4.5 

4.5 

3.5 

3.5 

4 

2.5 

1.75 

1 

3 

5 

2 

2.25 

1 

3 

Table  B.6:  Wavelet:  Real  Morlet,  SNR:3  dB 


s 

Real 

Imag. 

Angle 

Abs. 

1 

3 

1.75 

1 

3.75 

2 

4.5 

4 

2.75 

3.75 

3 

4.5 

4.5 

3 

4 

4 

2.75 

2.75 

1 

2 

5 

2.75 

2.75 

1 

1.5 

5  —  sharply  obvious  with  sharp  boundaries. 

4  —  obvious. 

3  —  fair. 

2  —  distinguishable  from  background  noise. 

1  —  not  distinguishable  from  background  noise. 
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Table  B.7:  Wavelet:  Complex  Daub-3,  SNR:3  dB 


s 

Real 

Imag. 

Angle 

Abs. 

1 

2.8 

3.6 

2.2 

1 

2 

2.4 

3.6 

4.25 

2.8 

3 

3.4 

4.25 

1.4 

1 

4 

1.8 

1 

1 

1 

5 

1 

1 

1 

1 

Table  B.8:  Wavelet:  Real  Daub-3,  SNR:3  dB 


s 

Real 

Imag. 

Angle 

Abs. 

1 

3.8 

1.4 

2.4 

1.6 

2 

3.6 

3.2 

4.4 

3.4 

3 

2.6 

3.2 

2.2 

1.2 

4 

1 

1 

1 

1 

5 

1 

1 

1 

1 

5  —  sharply  obvious  with  sharp  boundaries. 

4  —  obvious. 

3  —  fair. 

2  —  distinguishable  from  background  noise. 

1  —  not  distinguishable  from  background  noise. 
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APPENDIX  C.  WAVELET  SURFACES  OF 

OTHER  SIGNALS 

This  appendix  includes  the  wavelet  surfaces  obtained  by  processing  noise  only, 
ASK,  PSK  and  FSK  signals.  Figures  40  and  41  show  no  identifiable  diamond  structure 
pattern  for  the  noise  only  case.  Each  pair  of  successive  figures  (Figure  42  and  43, 
Figure  44  and  45,  and  Figure  46  and  47)  show  a  diamond  structure  residing  only  in 
one  scale  . 
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Figure  40.  Wavelet  surfaces  obtained  from  noise  only,  scale  1,  2  and  3. 
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Figure  41.  Wavelet  surfaces  obtained  from  noise  only,  scale  4  and  5. 
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Figure  42.  Wavelet  surfaces  obtained  from  ASK  signal,  scale  1,  2,  and  3. 
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Figure  43.  Wavelet  surfaces  obtained  from  ASK  signal,  scale  4  and  5. 
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Figure  44.  Wavelet  surfaces  obtained  from  FSK  signal,  scale  1,  2  and  3. 
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Figure  45.  Wavelet  surfaces  obtained  from  FSK  signal,  scale  4  and  5. 
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Figure  46.  Wavelet  surfaces  obtained  from  PSK  signal,  scale  1,  2  and  3. 
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Figure  47.  Wavelet  surfaces  obtained  from  PSK  signal,  scale  4  and  5. 
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